When not to avoid inbreeding: a gene’s eye view perspective

Load packages

# for tidy style coding and plotting

library(tidyverse) 
library(vroom) # to read lots of csv files at once

# more table options

library(pander) # for tables
library(kableExtra) # for scrolling tables
library(data.table) # for efficient handling of large dataframes

# making ggplot more powerful

library(MetBrewer) # for colour palettes based upon artwork housed at the MET
library(MoMAColors) # for colour palettes based upon artwork housed at MoMA
library(wesanderson) # for colour palettes based on wes anderson movies
library(rcartocolor) # for nice sequential colour schemes
library(PNWColors) # for colour palettes 
library(tidybayes) # for plotting distributions
library(stickylabeller) # labelling facets with strings in ggplot
library(geomtextpath) # for curved plot annotations
library(ggtext) # for markdown syntax in plot labels
library(patchwork) # for patching plots together
library(ggnewscale) # to reset scales in plots, allowing multiple fill arguments in ggplot

# for computation speed checks

library(profvis) # breakdown of complex functions
library(bench) # individual functions

PLotting results from the inclusive fitness model

We find the inbreeding tolerance thresholds with strict polygyny are

\[\begin{align} \delta_{\mathrm{f},i}^* = \frac{rc_{\mathrm{m},i}F_\mathrm{m} + rc_{\mathrm{m},i}}{c_{\mathrm{f},i}F_\mathrm{f} + rc_{\mathrm{m},i}F_\mathrm{m} + rc_{\mathrm{m},i} + c_{\mathrm{f},i}} \label{female_new_threshold} \end{align}\]

\[\begin{align} \delta_{\mathrm{m},i}^* = \frac{c_{\mathrm{m},i} + c_{\mathrm{m},i}F_\mathrm{m}}{c_{\mathrm{m},i} + c_{\mathrm{m},i}F_\mathrm{m} + rc_{\mathrm{f},i} + rc_{\mathrm{f},i}F_\mathrm{f}}. \label{male_new_threshold} \end{align}\]

We use these equations to produce Figure 1

Prepare for plotting

resolution <- 200

analytical_results <-
    expand_grid(
    r = seq(0, 1, length = resolution),
    c_m = c(0, 0.5, 1),
    c_f = c(0, 0.5, 1),
    F_m = c(0, 1),
    F_f = c(0, 1),
    D = seq(0, 1, length = resolution)) %>%  # D represents inbreeding depression) %>% 
  mutate(female_inbreeding_fitness = c_f*((1+F_f)/2)*(1-D) + (r*c_m*((1+F_m)/2)*(1-D)),
         male_inbreeding_fitness = c_m*((1+F_m)/2)*(1-D) + r*c_f*((1+F_f)/2)*(1-D),
         female_outbreeding_fitness = c_f*((1+F_f)/2),
         male_outbreeding_fitness = r*c_f*((1+F_f)/2),
         female_fitness_contrast = female_inbreeding_fitness - female_outbreeding_fitness, # this is close to a selection coefficient
         male_fitness_contrast = male_inbreeding_fitness - male_outbreeding_fitness,
         qual_f_contrast = if_else(female_fitness_contrast > 0, 1, -1),
         qual_m_contrast = if_else(male_fitness_contrast > 0, 1, -1))

# Extract results for females  
  
Within_organism_conflict_females <-
  analytical_results %>%
  select(1:6, female_fitness_contrast, qual_f_contrast) %>%
  mutate(Case = case_when(c_m == 1 & c_f == 1 & F_m == F_f ~ "Autosome",
                          c_m == 0.5 & c_f == 1 & F_m == 1 ~ "X",
                          c_m == 1 & c_f == 0.5 & F_f == 1 ~ "Z",
                          c_m == 0 & c_f == 1 & F_m == 0 & F_f == 1 ~ "W",
                          c_m == 0 & c_f == 1 & F_m == 1 & F_f == 1 ~ "Cytoplasmic (m)",
                          c_m == 1 & c_f == 0 & F_m == 1 & F_f == 1 ~ "Cytoplasmic (p)"),
         Inbreeding_speed = case_when(Case == "X" & F_f == 0 |
                                        Case == "Z" & F_m == 0 |
                                        Case == "Autosome" & F_f == 0 ~ "Slow",
                                      .default = "Fast")) %>% 
  filter(!is.na(Case),
         Inbreeding_speed == "Fast") %>%
  select(-c(c_f, c_m, F_f, F_m, Inbreeding_speed)) %>% 
  pivot_wider(names_from = Case, values_from = qual_f_contrast, id_cols = c(r, D)) %>% 
  mutate(conflict = case_when(Z < 0 ~ 1,
                              Autosome < 0 ~ 2,
                              X < 0 ~ 3,
                              W < 0 & `Cytoplasmic (m)` < 0 ~ 4))

# Extract results for males 

Within_organism_conflict_males <-
  analytical_results %>%
  select(1:6, male_fitness_contrast, qual_m_contrast) %>%
  mutate(Case = case_when(c_m == 1 & c_f == 1 & F_m == F_f ~ "Autosome",
                          c_m == 0.5 & c_f == 1 & F_m == 1 ~ "X",
                          c_m == 1 & c_f == 0.5 & F_f == 1 ~ "Z",
                          c_m == 1 & c_f == 0 & F_m == 1 & F_f == 0 ~ "Y",
                          c_m == 0 & c_f == 1 & F_m == 1 & F_f == 1 ~ "Cytoplasmic (m)",
                          c_m == 1 & c_f == 0 & F_m == 1 & F_f == 1 ~ "Cytoplasmic (p)"),
         Inbreeding_speed = case_when(Case == "X" & F_f == 0 |
                                        Case == "Z" & F_m == 0 |
                                        Case == "Autosome" & F_f == 0 ~ "Slow",
                                      .default = "Fast")) %>% 
  filter(!is.na(Case),
         Inbreeding_speed == "Fast") %>%
  select(-c(c_f, c_m, F_f, F_m, Inbreeding_speed)) %>% 
  pivot_wider(names_from = Case, values_from = qual_m_contrast, id_cols = c(r, D)) %>% 
  mutate(conflict = case_when(Z < 0 ~ 1,
                              Autosome < 0 ~ 2,
                              X < 0 ~ 3,
                              `Cytoplasmic (m)` < 0 ~ 4))

# calculate the inbreeding thresholds for the plot

figure_1_thresholds <- 
  expand_grid(r = seq(from = 0, to = 1, by = 0.01),
              c_m = c(0, 0.5, 1),
              c_f = c(0, 0.5, 1),
              F_m = c(0, 1),
              F_f = c(0, 1),
              Sex = c("Female", "Male")) %>% 
  mutate(D = case_when(Sex == "Female" ~ (r*c_m*F_m + r*c_m) / 
                         (r*c_m*F_m + r*c_m + c_f*F_f + c_f),
                       Sex == "Male" ~ (c_m + c_m*F_m) / 
                         (c_m + c_m*F_m + r*c_f + r*c_f*F_f))) %>% 
  mutate(Case = case_when(c_m == 1 & c_f == 1 & F_m == F_f ~ "Autosome",
                          c_m == 0.5 & c_f == 1 & F_m == 1 ~ "X",
                          c_m == 1 & c_f == 0.5 & F_f == 1 ~ "Z",
                          c_m == 0 & c_f == 1 & F_m == 0 & F_f == 1 & Sex == "Female" ~ "W & Cytoplasmic (m)",
                          c_m == 1 & c_f == 0 & F_m == 1 & F_f == 0 & Sex == "Male"~ "Y & Cytoplasmic (p)",
                          c_m == 0 & c_f == 1 & F_m == 1 & F_f == 1 & Sex == "Male" ~ "Cytoplasmic (m)",
                          c_m == 1 & c_f == 0 & F_m == 1 & F_f == 1 & Sex == "Female" ~ "Cytoplasmic (p)"),
         Inbreeding_speed = case_when(Case == "X" & F_f == 0 |
                                        Case == "Z" & F_m == 0 |
                                        Case == "Autosome" & F_m == 0 ~ "Slow",
                                      .default = "Fast")) %>% 
  filter(!is.na(Case),
         Inbreeding_speed == "Fast")

threshold_1 <- 
  figure_1_thresholds %>% 
  # remove labels that don't plot well
  filter(Case != "Y & Cytoplasmic (p)",
         Case != "Cytoplasmic (p)")

threshold_2 <- 
  figure_1_thresholds %>% 
  # remove labels that don't plot well
  filter(Case == "Y & Cytoplasmic (p)" |
         Case == "Cytoplasmic (p)")

Let’s plot the conflict zones between chromosomes, within the sexes. These form panels A and B of Figure 1.

greens <- c("#4c9b82", "#6cc08b", "#97e196", "#d3f2a3") 

#d3f2a3,#97e196,#6cc08b,#4c9b82,#217a79,#105965,#074050

fig_1a <-
  Within_organism_conflict_females %>% 
  ggplot(aes(x = r, y = D)) +
  geom_blank() +
  geom_tile(aes(fill = conflict)) +
  scale_fill_gradientn(colours = greens) +
  guides(fill="none") +
  geom_textline(data = threshold_1 %>% filter(Sex == "Female"), 
                aes(x = r, y = D, group = Case, label = Case), size = 4, vjust = -0.25) +
  geom_textline(data = threshold_2 %>% filter(Sex == "Female"), 
                aes(x = r, y = D, group = Case, label = Case), size = 4, vjust = 1.2) +
  coord_fixed() +
  labs(y = ~delta~'(inbreeding depression)', 
       x = '_r_ (Individual-level relatedness)',
       title = "Female within-organism conflict") +
  scale_x_continuous(expand = c(0, 0), 
                     breaks = c(0, 0.25, 0.5, 0.75, 1)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.title.x = element_markdown(size =14),
        axis.title.y = element_text(size=14),
        axis.text = element_text(size = 12),
        plot.title = element_text(hjust = 0.5, size = 14))

fig_1b <-
  Within_organism_conflict_males %>% 
  ggplot(aes(x = r, y = D)) +
  geom_blank() +
  geom_tile(aes(fill = conflict)) +
  scale_fill_gradientn(colours = greens) +
  guides(fill="none") +
  geom_textline(data = threshold_1 %>% filter(Sex == "Male"), 
                aes(x = r, y = D, group = Case, label = Case), size = 4, vjust = -0.25) +
  geom_textline(data = threshold_2 %>% filter(Sex == "Male"), 
                aes(x = r, y = D, group = Case, label = Case), size = 4, vjust = 1.2) +
  coord_fixed() +
  labs(y = ~delta~'(inbreeding depression)', 
       x = '_r_ (Individual-level relatedness)',
       title = "Male within-organism conflict") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.25, 0.5, 0.75, 1)) + 
  scale_y_continuous(expand = c(0, 0)) + # labels = c(0, 25, 50, 75, 90)) +
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.title.x = element_markdown(size =14),
        axis.title.y = element_text(size=14),
        axis.text = element_text(size = 12),
        plot.title = element_text(hjust = 0.5, size = 14))

top <-
  fig_1a + fig_1b & plot_layout(guides = 'collect', axes = "collect")

top

Figure 1A-B. Within-organism conflict over the expression of inbreeding. In panels and , contours show the the severity of inbreeding depression below which inbreeding increases the inclusive fitness of genes in a particular genomic region. Results are shown for the fast inbreeding scenario (see main text). As the fill gets darker, fewer parts of the genome have higher fitness when inbreeding is preferred over outcrossing. Cytoplasmic (m) and (p) denote cytoplasmic genes that have strict maternal and paternal inheritance, respectively.

Make Figure 1C-D

In the bottom two panels of Figure 1, we show offspring viabilities are the inbreeding tolerance threshold for alleles in each genomic region.

fitness_loss_f <-
  figure_1_thresholds %>% 
  filter(r == 0.5) %>% 
  mutate(Case = as.factor(Case),
         Sex = case_when(Sex == "Female" ~ "delta[f]~when~r==0.5",
                         .default = "delta[m]~when~r==0.5")) %>% 
  filter(Sex == "delta[f]~when~r==0.5") %>% 
  mutate(chromosome = fct_relevel(Case, c("Autosome", "X", "Z", "W & Cytoplasmic (m)",
                                          "Cytoplasmic (p)"))) %>% 
  ggplot(aes(x = reorder(Case, D), y = 1 - D)) +
  geom_col(fill = "#d0e2af", alpha = 0.75) +
  facet_wrap(~Sex, labeller = label_parsed) +
  scale_y_continuous(expand = c(0.0, 0.0), limits = c(0, 1.01), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + 
  labs(x = "Genomic location",
       y = "Offspring viability relative\nto outbreeding population") +
  theme_bw() +
  theme(strip.background = element_rect(colour = "black", fill = "#d0e2af"),
        strip.text = element_text(size = 14),
        axis.title = element_text(size = 14),
        axis.text.y = element_text(size = 12))

fitness_loss_m <-
  figure_1_thresholds %>% 
  filter(r == 0.5) %>% 
  mutate(Case = as.factor(Case),
         Sex = case_when(Sex == "Female" ~ "delta[f]~when~r==0.5",
                         .default = "delta[m]~when~r==0.5")) %>% 
  filter(Sex == "delta[m]~when~r==0.5") %>% 
  mutate(chromosome = fct_relevel(Case, c("Autosome", "X", "Z", "Y & Cytoplasmic (p)",
                                          "Cytoplasmic (m)"))) %>% 
  ggplot(aes(x = reorder(Case, D), y = 1 - D)) +
  geom_col(fill = "#d0e2af", alpha = 0.75) +
  facet_wrap(~Sex, labeller = label_parsed) +
  scale_y_continuous(expand = c(0.0, 0.0), limits = c(0, 1.01), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + 
  labs(x = "Genomic location",
       y = "Offspring viability relative\nto outbreeding population") +
  theme_bw() +
  theme(strip.background = element_rect(colour = "black", fill = "#d0e2af"),
        strip.text = element_text(size = 14),
        axis.title = element_text(size = 14),
        axis.text.y = element_text(size = 12))
  

    

bottom <- 
  fitness_loss_f + fitness_loss_m +
  plot_layout(axis_titles = "collect")

bottom

Make Figure 1

top / bottom + plot_annotation(tag_levels = 'A')

Figure 1. Within-organism conflict over the expression of inbreeding. In panels and , contours show the the severity of inbreeding depression below which inbreeding increases the inclusive fitness of genes in a particular genomic region. Results are shown for the fast inbreeding scenario. As the fill gets darker, fewer parts of the genome have higher fitness when inbreeding is preferred over outcrossing. and show the lowest offspring viabilities that permit the invasion of inbreeding alleles when mating is between full-siblings. Cytoplasmic (m) and (p) denote cytoplasmic genes that have strict maternal and paternal inheritance, respectively.

Make Figure SX

# calculate every combination

Intragenomic_all_combos <-
  analytical_results %>%
  select(1:5, female_fitness_contrast, male_fitness_contrast) %>%
  mutate(Case = case_when(a_m == 0 & a_f == 0.5 ~ "W",
                          a_m == 0 & a_f == 1 ~ "Cytoplasmic",
                          a_m == 0.5 & a_f == 1 ~ "X",
                          a_m == 0.5 & a_f == 0 ~ "Y",
                          a_m == 1 & a_f == 0.5 ~ "Z",
                          a_m == 1 & a_f == 1 ~ "Autosome")) %>% 
  filter(!is.na(Case) & c == 0) %>%
  select(-c(a_f, a_m)) %>% 
  pivot_wider(names_from = Case, 
              values_from = c(female_fitness_contrast, male_fitness_contrast)) %>% 
  mutate(`Autosome (female) vs Autosome (male)` = I_calculator(female_fitness_contrast_Autosome, male_fitness_contrast_Autosome),
         `Autosome (female) vs X (female)` = I_calculator(female_fitness_contrast_Autosome,
                                                          female_fitness_contrast_X),
         `Autosome (female) vs X (male)` = I_calculator(female_fitness_contrast_Autosome,
                                                        male_fitness_contrast_X),
         `Autosome (female) vs Y (male)` = I_calculator(female_fitness_contrast_Autosome,
                                                        male_fitness_contrast_Y),
         `Autosome (female) vs Z (female)` = I_calculator(female_fitness_contrast_Autosome,
                                                          female_fitness_contrast_Z),
         `Autosome (female) vs Z (male)` = I_calculator(female_fitness_contrast_Autosome,
                                                        male_fitness_contrast_Z),
         `Autosome (female) vs W or Cytoplasmic (female)` = I_calculator(female_fitness_contrast_Autosome,
                                                                         `female_fitness_contrast_W or Cytoplasmic`),
         `Autosome (female) vs Cytoplasmic (male)` = I_calculator(female_fitness_contrast_Autosome,
                                                                  `male_fitness_contrast_W or Cytoplasmic`),
         
         `X (female) vs Autosome (male)` = I_calculator(female_fitness_contrast_X, male_fitness_contrast_Autosome),
         `X (female) vs X (male)` = I_calculator(female_fitness_contrast_X, male_fitness_contrast_X),
         `X (female) vs Y (male)` = I_calculator(female_fitness_contrast_X, male_fitness_contrast_Y),
         `X (female) vs Z (female)` = I_calculator(female_fitness_contrast_X, female_fitness_contrast_Z),
         `X (female) vs Z (male)` = I_calculator(female_fitness_contrast_X, male_fitness_contrast_Z),
         `X (female) vs W or Cytoplasmic (female)` = I_calculator(female_fitness_contrast_X,
                                                                  `female_fitness_contrast_W or Cytoplasmic`),
         `X (female) vs Cytoplasmic (male)` = I_calculator(female_fitness_contrast_X,
                                                           `male_fitness_contrast_W or Cytoplasmic`),
         
         `Z (female) vs Autosome (male)` = I_calculator(female_fitness_contrast_Z, male_fitness_contrast_Autosome),
         `Z (female) vs X (male)` = I_calculator(female_fitness_contrast_Z, male_fitness_contrast_X),
         `Z (female) vs Y (male)` = I_calculator(female_fitness_contrast_Z, male_fitness_contrast_Y),
         `Z (female) vs Z (male)` = I_calculator(female_fitness_contrast_Z, male_fitness_contrast_Z),
         `Z (female) vs W or Cytoplasmic (female)` = I_calculator(female_fitness_contrast_Z, 
                                                                  `female_fitness_contrast_W or Cytoplasmic`),
         `Z (female) vs Cytoplasmic (male)` = I_calculator(female_fitness_contrast_Z, 
                                                           `male_fitness_contrast_W or Cytoplasmic`),
         
         `W or Cytoplasmic (female) vs Autosome (male)` = I_calculator(`female_fitness_contrast_W or Cytoplasmic`,
                                                                       male_fitness_contrast_Autosome),
         `W or Cytoplasmic (female) vs X (male)` = I_calculator(`female_fitness_contrast_W or Cytoplasmic`,
                                                                male_fitness_contrast_X),
         `W or Cytoplasmic (female) vs Y (male)` = I_calculator(`female_fitness_contrast_W or Cytoplasmic`,
                                                                male_fitness_contrast_Y),
         `W or Cytoplasmic (female) vs Z (male)` = I_calculator(`female_fitness_contrast_W or Cytoplasmic`,
                                                                male_fitness_contrast_Z),
         `W or Cytoplasmic (female) vs Cytoplasmic (male)` = I_calculator(`female_fitness_contrast_W or Cytoplasmic`,
                                                                          `male_fitness_contrast_W or Cytoplasmic`),
         
         `Autosome (male) vs X (male)` = I_calculator(male_fitness_contrast_Autosome, 
                                                      male_fitness_contrast_X),
         `Autosome (male) vs Y (male)` = I_calculator(male_fitness_contrast_Autosome, 
                                                      male_fitness_contrast_Y),
         `Autosome (male) vs Z (male)` = I_calculator(male_fitness_contrast_Autosome, 
                                                      male_fitness_contrast_Z),
         `Autosome (male) vs Cytoplasmic (male)` = I_calculator(male_fitness_contrast_Autosome, 
                                                                `male_fitness_contrast_W or Cytoplasmic`),
         
         `X (male) vs Y (male)` = I_calculator(male_fitness_contrast_X, male_fitness_contrast_Y),
         `X (male) vs Z (male)` = I_calculator(male_fitness_contrast_X, male_fitness_contrast_Z),
         `X (male) vs Cytoplasmic (male)` = I_calculator(male_fitness_contrast_X, 
                                                         `male_fitness_contrast_W or Cytoplasmic`),
         
         `Y (male) vs Z (male)` = I_calculator(male_fitness_contrast_Y, male_fitness_contrast_Z),
         `Y (male) vs Cytoplasmic (male)` = I_calculator(male_fitness_contrast_Y, 
                                                         `male_fitness_contrast_W or Cytoplasmic`),
         
         `Z (male) vs Cytoplasmic (male)` = I_calculator(male_fitness_contrast_Z, 
                                                         `male_fitness_contrast_W or Cytoplasmic`)
  ) %>%
  pivot_longer(cols = 14:49, names_to = "contrast", values_to = "Evolutionary_conflict") %>% 
  mutate(relationship = case_when(
    # beneficial in both contexts
    contrast == "Autosome (female) vs Autosome (male)" 
    & female_fitness_contrast_Autosome > 0 & male_fitness_contrast_Autosome > 0 | 
      contrast == "Autosome (female) vs X (female)" &
      female_fitness_contrast_Autosome > 0 & female_fitness_contrast_X > 0 | 
      contrast == "Autosome (female) vs X (male)" &
      female_fitness_contrast_Autosome > 0 & male_fitness_contrast_X > 0 |
      contrast == "Autosome (female) vs Y (male)" &
      female_fitness_contrast_Autosome > 0 & male_fitness_contrast_Y > 0 |
      contrast == "Autosome (female) vs Z (female)" &
      female_fitness_contrast_Autosome > 0 & female_fitness_contrast_Z > 0 |
      contrast == "Autosome (female) vs Z (male)" &
      female_fitness_contrast_Autosome > 0 & male_fitness_contrast_Z > 0 | 
      contrast == "Autosome (female) vs W or Cytoplasmic (female)" &
      female_fitness_contrast_Autosome > 0 & `female_fitness_contrast_W or Cytoplasmic` > 0 | 
      contrast == "Autosome (female) vs Cytoplasmic (male)"  &
      female_fitness_contrast_Autosome > 0 & `male_fitness_contrast_W or Cytoplasmic` > 0 |
      
      contrast == "X (female) vs Autosome (male)" & 
      female_fitness_contrast_X > 0 & male_fitness_contrast_Autosome > 0 | 
      contrast == "X (female) vs X (male)" & female_fitness_contrast_X > 0 & male_fitness_contrast_X > 0 |
      contrast == "X (female) vs Y (male)" & female_fitness_contrast_X > 0 & male_fitness_contrast_Y > 0 |
      contrast == "X (female) vs Z (female)" & female_fitness_contrast_X > 0  & female_fitness_contrast_Z > 0|
      contrast == "X (female) vs Z (male)" & female_fitness_contrast_X > 0 & male_fitness_contrast_Z > 0 | 
      contrast == "X (female) vs W or Cytoplasmic (female)" &
      female_fitness_contrast_X > 0 & `female_fitness_contrast_W or Cytoplasmic` > 0 |
      contrast == "X (female) vs Cytoplasmic (male)" &
      female_fitness_contrast_X > 0 & `male_fitness_contrast_W or Cytoplasmic` > 0 |
      
      contrast == "Z (female) vs Autosome (male)" &
      female_fitness_contrast_Z > 0 & male_fitness_contrast_Autosome > 0 | 
      contrast == "Z (female) vs X (male)" & female_fitness_contrast_Z > 0 & male_fitness_contrast_X > 0 |
      contrast == "Z (female) vs Y (male)" & female_fitness_contrast_Z > 0 & male_fitness_contrast_Y > 0 |
      contrast == "Z (female) vs Z (male)" & female_fitness_contrast_Z > 0 & male_fitness_contrast_Z > 0 | 
      contrast == "Z (female) vs W or Cytoplasmic (female)" &
      female_fitness_contrast_Z > 0 & `female_fitness_contrast_W or Cytoplasmic` > 0 | 
      contrast == "Z (female) vs Cytoplasmic (male)" &
      female_fitness_contrast_Z > 0 & `male_fitness_contrast_W or Cytoplasmic` > 0 | 
      
      contrast == "W or Cytoplasmic (female) vs Autosome (male)" &
      `female_fitness_contrast_W or Cytoplasmic` > 0 & male_fitness_contrast_Autosome > 0 | 
      contrast == "W or Cytoplasmic (female) vs X (male)" &
      `female_fitness_contrast_W or Cytoplasmic` > 0 & male_fitness_contrast_X > 0 |
      contrast == "W or Cytoplasmic (female) vs Y (male)" &
      `female_fitness_contrast_W or Cytoplasmic` > 0 & male_fitness_contrast_Y > 0 |
      contrast == "W or Cytoplasmic (female) vs Z (male)" &
      `female_fitness_contrast_W or Cytoplasmic` > 0 & male_fitness_contrast_Z > 0 | 
      contrast == "W or Cytoplasmic (female) vs Cytoplasmic (male)" &
      `female_fitness_contrast_W or Cytoplasmic` > 0 & `male_fitness_contrast_W or Cytoplasmic` > 0 | 
      
      contrast == "Autosome (male) vs X (male)" & male_fitness_contrast_Autosome > 0 & male_fitness_contrast_X > 0 | 
      contrast == "Autosome (male) vs Y (male)" & male_fitness_contrast_Autosome > 0 & male_fitness_contrast_Y > 0 |
      contrast == "Autosome (male) vs Z (male)" & male_fitness_contrast_Autosome > 0 & male_fitness_contrast_Z > 0 |
      contrast == "Autosome (male) vs Cytoplasmic (male)" &
      male_fitness_contrast_Autosome > 0 & `male_fitness_contrast_W or Cytoplasmic` > 0 |
      
      contrast == "X (male) vs Y (male)" & male_fitness_contrast_X > 0 & male_fitness_contrast_Y > 0 |
      contrast == "X (male) vs Z (male)" & male_fitness_contrast_X > 0 & male_fitness_contrast_Z > 0 |
      contrast == "X (male) vs Cytoplasmic (male)" &
      male_fitness_contrast_X > 0 & `male_fitness_contrast_W or Cytoplasmic` > 0 |
      
      contrast == "Y (male) vs Z (male)" & male_fitness_contrast_Y > 0 & male_fitness_contrast_Z > 0 |
      contrast == "Y (male) vs Cytoplasmic (male)" &
      male_fitness_contrast_Y > 0 & `male_fitness_contrast_W or Cytoplasmic` > 0 |
      
      contrast == "Z (male) vs Cytoplasmic (male)" &
      male_fitness_contrast_Z > 0 & `male_fitness_contrast_W or Cytoplasmic` > 0
      
      ~ "Inbreeding favoured in both contexts",
    
    # deleterious in both contexts
    
       contrast == "Autosome (female) vs Autosome (male)" 
    & female_fitness_contrast_Autosome < 0 & male_fitness_contrast_Autosome < 0 | 
      contrast == "Autosome (female) vs X (female)" &
      female_fitness_contrast_Autosome < 0 & female_fitness_contrast_X < 0 | 
      contrast == "Autosome (female) vs X (male)" &
      female_fitness_contrast_Autosome < 0 & male_fitness_contrast_X < 0 |
      contrast == "Autosome (female) vs Y (male)" &
      female_fitness_contrast_Autosome < 0 & male_fitness_contrast_Y < 0 |
      contrast == "Autosome (female) vs Z (female)" &
      female_fitness_contrast_Autosome < 0 & female_fitness_contrast_Z < 0 |
      contrast == "Autosome (female) vs Z (male)" &
      female_fitness_contrast_Autosome < 0 & male_fitness_contrast_Z < 0 | 
      contrast == "Autosome (female) vs W or Cytoplasmic (female)" &
      female_fitness_contrast_Autosome < 0 & `female_fitness_contrast_W or Cytoplasmic` < 0 | 
      contrast == "Autosome (female) vs Cytoplasmic (male)"  &
      female_fitness_contrast_Autosome < 0 & `male_fitness_contrast_W or Cytoplasmic` < 0 |
      
      contrast == "X (female) vs Autosome (male)" & 
      female_fitness_contrast_X < 0 & male_fitness_contrast_Autosome < 0 | 
      contrast == "X (female) vs X (male)" & female_fitness_contrast_X < 0 & male_fitness_contrast_X < 0 |
      contrast == "X (female) vs Y (male)" & female_fitness_contrast_X < 0 & male_fitness_contrast_Y < 0 |
      contrast == "X (female) vs Z (female)" & female_fitness_contrast_X < 0  & female_fitness_contrast_Z < 0|
      contrast == "X (female) vs Z (male)" & female_fitness_contrast_X < 0 & male_fitness_contrast_Z < 0 | 
      contrast == "X (female) vs W or Cytoplasmic (female)" &
      female_fitness_contrast_X < 0 & `female_fitness_contrast_W or Cytoplasmic` < 0 |
      contrast == "X (female) vs Cytoplasmic (male)" &
      female_fitness_contrast_X < 0 & `male_fitness_contrast_W or Cytoplasmic` < 0 |
      
      contrast == "Z (female) vs Autosome (male)" &
      female_fitness_contrast_Z < 0 & male_fitness_contrast_Autosome < 0 | 
      contrast == "Z (female) vs X (male)" & female_fitness_contrast_Z < 0 & male_fitness_contrast_X < 0 |
      contrast == "Z (female) vs Y (male)" & female_fitness_contrast_Z < 0 & male_fitness_contrast_Y < 0 |
      contrast == "Z (female) vs Z (male)" & female_fitness_contrast_Z < 0 & male_fitness_contrast_Z < 0 | 
      contrast == "Z (female) vs W or Cytoplasmic (female)" &
      female_fitness_contrast_Z < 0 & `female_fitness_contrast_W or Cytoplasmic` < 0 | 
      contrast == "Z (female) vs Cytoplasmic (male)" &
      female_fitness_contrast_Z < 0 & `male_fitness_contrast_W or Cytoplasmic` < 0 | 
      
      contrast == "W or Cytoplasmic (female) vs Autosome (male)" &
      `female_fitness_contrast_W or Cytoplasmic` < 0 & male_fitness_contrast_Autosome < 0 | 
      contrast == "W or Cytoplasmic (female) vs X (male)" &
      `female_fitness_contrast_W or Cytoplasmic` < 0 & male_fitness_contrast_X < 0 |
      contrast == "W or Cytoplasmic (female) vs Y (male)" &
      `female_fitness_contrast_W or Cytoplasmic` < 0 & male_fitness_contrast_Y < 0 |
      contrast == "W or Cytoplasmic (female) vs Z (male)" &
      `female_fitness_contrast_W or Cytoplasmic` < 0 & male_fitness_contrast_Z < 0 | 
      contrast == "W or Cytoplasmic (female) vs Cytoplasmic (male)" &
      `female_fitness_contrast_W or Cytoplasmic` < 0 & `male_fitness_contrast_W or Cytoplasmic` < 0 | 
      
      contrast == "Autosome (male) vs X (male)" & male_fitness_contrast_Autosome < 0 & male_fitness_contrast_X < 0 | 
      contrast == "Autosome (male) vs Y (male)" & male_fitness_contrast_Autosome < 0 & male_fitness_contrast_Y < 0 |
      contrast == "Autosome (male) vs Z (male)" & male_fitness_contrast_Autosome < 0 & male_fitness_contrast_Z < 0 |
      contrast == "Autosome (male) vs Cytoplasmic (male)" &
      male_fitness_contrast_Autosome < 0 & `male_fitness_contrast_W or Cytoplasmic` < 0 |
      
      contrast == "X (male) vs Y (male)" & male_fitness_contrast_X < 0 & male_fitness_contrast_Y < 0 |
      contrast == "X (male) vs Z (male)" & male_fitness_contrast_X < 0 & male_fitness_contrast_Z < 0 |
      contrast == "X (male) vs Cytoplasmic (male)" &
      male_fitness_contrast_X < 0 & `male_fitness_contrast_W or Cytoplasmic` < 0 |
      
      contrast == "Y (male) vs Z (male)" & male_fitness_contrast_Y < 0 & male_fitness_contrast_Z < 0 |
      contrast == "Y (male) vs Cytoplasmic (male)" &
      male_fitness_contrast_Y < 0 & `male_fitness_contrast_W or Cytoplasmic` < 0 |
      
      contrast == "Z (male) vs Cytoplasmic (male)" &
      male_fitness_contrast_Z < 0 & `male_fitness_contrast_W or Cytoplasmic` < 0
    
    ~ "Inbreeding deleterious in both contexts",
    
    .default = "Intragenomic conflict")) %>% 
  mutate(Evolutionary_conflict = if_else(is.na(Evolutionary_conflict), 0, Evolutionary_conflict),
         contrast = fct_relevel(contrast, 
                                "W or Cytoplasmic (female) vs Cytoplasmic (male)",
                                "X (female) vs Cytoplasmic (male)", 
                                "X (female) vs W or Cytoplasmic (female)",
                                "Autosome (female) vs Cytoplasmic (male)",
                                "Autosome (female) vs W or Cytoplasmic (female)",
                                "Z (female) vs Cytoplasmic (male)",
                                "Z (female) vs W or Cytoplasmic (female)",
                                "X (male) vs Cytoplasmic (male)",
                                "W or Cytoplasmic (female) vs X (male)",
                                "Autosome (male) vs Cytoplasmic (male)",
                                "W or Cytoplasmic (female) vs Autosome (male)",
                                "Z (male) vs Cytoplasmic (male)",
                                "W or Cytoplasmic (female) vs Z (male)",
                                "Y (male) vs Cytoplasmic (male)",
                                "W or Cytoplasmic (female) vs Y (male)",
                                "Autosome (female) vs X (female)",
                                "X (female) vs Z (female)",
                                "Autosome (female) vs Z (female)",
                                "X (female) vs X (male)",
                                "X (female) vs Autosome (male)",
                                "Autosome (female) vs X (male)",
                                "X (female) vs Z (male)",
                                "Autosome (male) vs X (male)",
                                "Z (female) vs X (male)",
                                "Autosome (female) vs Autosome (male)",
                                "Autosome (female) vs Z (male)",
                                "Z (female) vs Autosome (male)",
                                "Z (female) vs Z (male)",
                                "X (male) vs Z (male)",
                                "Autosome (male) vs Z (male)",
                                "X (female) vs Y (male)",
                                "Autosome (female) vs Y (male)",
                                "Z (female) vs Y (male)",
                                "X (male) vs Y (male)",
                                "Autosome (male) vs Y (male)",
                                "Y (male) vs Z (male)"
                                ))

make_genomic_conflict_plot_2 <- 
  function(data, colour_pal){
    data %>% 
      ggplot(aes(x = r, y = D)) +
      geom_blank() +
      geom_tile(data = data %>% filter(relationship == "Intragenomic conflict"),
                aes(fill = Evolutionary_conflict*-1)) + 
      #geom_tile(data = data,
      #         aes(fill = Evolutionary_conflict)) + 
      scale_fill_gradientn(colours = colour_pal, limits = c(0, 0.6), #na.value = "white",
                           breaks = c(0, 0.6),
                           labels = c("Weaker conflict", "Stronger conflict")) +
      labs(fill = "Evolutionary conflict") +
      new_scale_fill() +
      geom_tile(data = data %>% filter(relationship != "Intragenomic conflict"),
                aes(fill = relationship), alpha = 0.75) +
      scale_fill_manual(values = c("#fbe6c5", "#d2fbd4"), 
                        labels = c("Inbreeding deleterious\nin both contexts", 
                                   "Inbreeding favoured\nin both contexts")) +
      stat_contour(aes(z = Evolutionary_conflict*-1), colour = "black", binwidth = 25,
                   breaks = c(0, 0.2, 0.4, 0.6)) +
      facet_wrap(~contrast, nrow = 6,
                 scales = "free", strip.position = c("top"),
                 labeller = label_glue('{`contrast`}')) +
      labs(x = "Individual-level relatedness coefficient",
           y = "Inbreeding depression",
           fill = "Evolutionary concordance") +
      scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.25, 0.5, 0.75)) + 
      scale_y_continuous(expand = c(0, 0)) + # labels = c(0, 25, 50, 75, 90)) +
      theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
            panel.grid.minor = element_blank(),
            strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
            strip.text = element_text(size = 6),
            axis.title.x = element_markdown(),
            axis.text = element_text(size = 8),
            #plot.title = element_text(hjust = 0.5, size = 8),
            legend.position = "none")
  }

make_genomic_conflict_plot_2(Intragenomic_all_combos, oranges)

Figure SX. all conflicts that can occur over inbreeding.

How fast do inbreeding families become homozygous?

Run a quick set of simulations to work this out, where female-male pairs of all possible genotypes at the inbreeding locus begin new families. We simulate the autosomal, X and Z cases, as these are where the number of I alleles carried by individuals matters.

p_space <-
  bind_rows(
    expand_grid(Dominance = c("dominant", "additive", "recessive"),
                starting_female = c(0, 1, 2),
                starting_male = c(0, 1, 2),
                expressed_sex = c("females", "males")) %>%  
      mutate(inheritance_scheme = "offspring_genotypes_autosome"),
    
    expand_grid(Dominance = c("dominant", "additive", "recessive"),
                starting_female = c(0, 1, 2),
                starting_male = c(0, 1, 2),
                expressed_sex = c("females", "males")) %>%  
      mutate(inheritance_scheme = "offspring_genotypes_X"),
    
    expand_grid(Dominance = c("dominant", "additive", "recessive"),
                starting_female = c(0, 1, 2),
                starting_male = c(0, 1, 2),
                expressed_sex = c("females", "males")) %>%  
      mutate(inheritance_scheme = "offspring_genotypes_Z"))  

family_homozygosity_sim <- function(row, # which parameter
                                    parameters, # the parameter space
                                    gens){  # max number of generations to run for
  
  expressed_sex <- parameters$expressed_sex[row]
  dominance <- parameters$Dominance[row]
  starting_female <- parameters$starting_female[row]
  starting_male <- parameters$starting_male[row]
  homozygous_genotype <- parameters$homozygous_genotype[row]
  heterozygous_genotype <- parameters$heterozygous_genotype[row]
  hemizygous_genotype <- parameters$hemizygous_genotype[row]
  inheritance_scheme <- eval(parse(text = parameters$inheritance_scheme[row])) 
  
  # start with individual pair
  mating_partners <- 
    data.table(Female_genotype = starting_female,
               Male_genotype = starting_male,
               mating_freq = 1)
  
  # set gen_counter to 1
  gen_counter <- 1
  
  tracking_table <- tibble(Female_homozygote_freq = numeric(), 
                               Female_heterozygote_freq = numeric(),
                               Female_hemizygote_freq = numeric(),
                               Male_homozygote_freq = numeric(), 
                               Male_heterozygote_freq = numeric(),
                               Male_hemizygote_freq = numeric(),
                               Generation = numeric()) # filled in as sim progresses
  
  mating_partners <- 
    data.table(Female_genotype = starting_female,
               Male_genotype = starting_male,
               mating_freq = 1)
  
  
  while(gen_counter <= gens){
    
    if(dominance == "dominant" & expressed_sex == "females"){
      offspring <- mating_partners[inheritance_scheme, 
                                   on = .(Female_genotype = Female_genotype,
                                          Male_genotype = Male_genotype), 
                                   nomatch = NULL, allow.cartesian  = TRUE
      ][, zygote_freq := zygote_freq * mating_freq
      ][Genotype == homozygous_genotype & Sex == 1 | # outbreeding females leave family
          Genotype == heterozygous_genotype & Sex == 1 | 
          Genotype == hemizygous_genotype & Sex == 1 |
          Sex == 0] 
    }
    
    if(dominance == "dominant" & expressed_sex == "males"){
      offspring <- mating_partners[inheritance_scheme, 
                                   on = .(Female_genotype = Female_genotype,
                                          Male_genotype = Male_genotype), 
                                   nomatch = NULL, allow.cartesian  = TRUE
      ][, zygote_freq := zygote_freq * mating_freq
      ][Genotype == homozygous_genotype & Sex == 0 | # outbreeding males leave family
          Genotype == heterozygous_genotype & Sex == 0 | 
          Genotype == hemizygous_genotype & Sex == 0 |
          Sex == 1] 
    }
    
    if(dominance == "additive" & expressed_sex == "females"){
      offspring <- mating_partners[inheritance_scheme, 
                                   on = .(Female_genotype = Female_genotype,
                                          Male_genotype = Male_genotype), 
                                   nomatch = NULL, allow.cartesian  = TRUE
      ][, zygote_freq := zygote_freq * mating_freq
      ][Genotype == heterozygous_genotype & Sex == 1, zygote_freq := zygote_freq * 0.5 # half the heterozygous females don't express I allele
      ][Genotype == homozygous_genotype & Sex == 1 | # outbreeding females leave family
          Genotype == heterozygous_genotype & Sex == 1 | 
          Genotype == hemizygous_genotype & Sex == 1 |
          Sex == 0] 
    }  
    
    if(dominance == "additive" & expressed_sex == "males"){
      offspring <- mating_partners[inheritance_scheme, 
                                   on = .(Female_genotype = Female_genotype,
                                          Male_genotype = Male_genotype), 
                                   nomatch = NULL, allow.cartesian  = TRUE
      ][, zygote_freq := zygote_freq * mating_freq
      ][Genotype == heterozygous_genotype & Sex == 0, zygote_freq := zygote_freq * 0.5 # half the heterozygous males don't express I allele
      ][Genotype == homozygous_genotype & Sex == 0 | # outbreeding males leave family
          Genotype == heterozygous_genotype & Sex == 0 | 
          Genotype == hemizygous_genotype & Sex == 0 |
          Sex == 1] 
    }    
    
    if(dominance == "recessive" & expressed_sex == "females"){
      offspring <- mating_partners[inheritance_scheme, 
                                   on = .(Female_genotype = Female_genotype,
                                          Male_genotype = Male_genotype), 
                                   nomatch = NULL, allow.cartesian  = TRUE
      ][, zygote_freq := zygote_freq * mating_freq
      ][Genotype == homozygous_genotype & Sex == 1 | # outbreeding females leave family
          Genotype == hemizygous_genotype & Sex == 1 |
          Sex == 0] 
    } 
    
    if(dominance == "recessive" & expressed_sex == "males"){
      offspring <- mating_partners[inheritance_scheme, 
                                   on = .(Female_genotype = Female_genotype,
                                          Male_genotype = Male_genotype), 
                                   nomatch = NULL, allow.cartesian  = TRUE
      ][, zygote_freq := zygote_freq * mating_freq
      ][Genotype == homozygous_genotype & Sex == 0 | # outbreeding females leave family
          Genotype == hemizygous_genotype & Sex == 0 |
          Sex == 1] 
    }   
    
    adult_females <- 
      offspring[Sex == 1, .(Female_genotype = Genotype, Female_freq = zygote_freq)
      ][, Female_freq := Female_freq / sum(Female_freq)
      ][, .(Female_freq = sum(Female_freq)), by = Female_genotype]
    
    if(nrow(adult_females[Female_genotype == homozygous_genotype]) > 0){ # fix this and the rest
      f_homozygote <- adult_females[Female_genotype == homozygous_genotype]$Female_freq
    } else f_homozygote <- 0
    
    if(nrow(adult_females[Female_genotype == heterozygous_genotype]) > 0){  
      f_heterozygote <- adult_females[Female_genotype == heterozygous_genotype]$Female_freq
    } else f_heterozygote <- 0
    
    if(nrow(adult_females[Female_genotype == hemizygous_genotype]) > 0){  
      f_hemizygote <- adult_females[Female_genotype == hemizygous_genotype]$Female_freq
    } else f_hemizygote <- 0
    
    adult_males <- 
      offspring[Sex == 0, .(Male_genotype = Genotype, Male_freq = zygote_freq)
      ][, Male_freq := Male_freq / sum(Male_freq)
      ][, .(Male_freq = sum(Male_freq)), by = Male_genotype]
    
    if(nrow(adult_males[Male_genotype == homozygous_genotype]) > 0){
      m_homozygote <- adult_males[Male_genotype == homozygous_genotype]$Male_freq
    } else m_homozygote <- 0
    
    if(nrow(adult_males[Male_genotype == heterozygous_genotype]) > 0){  
      m_heterozygote <- adult_males[Male_genotype == heterozygous_genotype]$Male_freq
    } else m_heterozygote <- 0
    
    if(nrow(adult_males[Male_genotype == hemizygous_genotype]) > 0){  
      m_hemizygote <- adult_males[Male_genotype == hemizygous_genotype]$Male_freq
    } else m_hemizygote <- 0
    
    tracking_table <- 
      rbindlist(list(tracking_table, 
                     list(f_homozygote,
                          f_heterozygote,
                          f_hemizygote,
                          m_homozygote,
                          m_heterozygote,
                          m_hemizygote,
                          gen_counter)), 
                ignore.attr=TRUE)
    
    mating_partners <- 
      CJ(Female_genotype = adult_females$Female_genotype, 
         Male_genotype = adult_males$Male_genotype)[
           adult_females, on = "Female_genotype"][
             adult_males, on = "Male_genotype"
           ][, mating_freq := Female_freq * Male_freq]
    
    gen_counter <- gen_counter + 1
    
  }
  tracking_table %>% 
    mutate(Dominance = dominance,
           starting_female = starting_female,
           starting_male = starting_male,
           expressed_sex = expressed_sex)
}

#family_homozygosity_sim(13, p_space, 10)

output <- map_dfr(1:nrow(p_space),
                      family_homozygosity_sim,
                      p_space,
                      20) 

plotting_data <-
  output %>% as_tibble() %>% 
  mutate(chromosome = case_when(str_detect(starting_female, "A") ~ "Autosome",
                                str_detect(starting_female, "X") ~ "X",
                                str_detect(starting_female, "Z") ~ "Z")) %>%
  unite("Mating pair", starting_female, starting_male, sep = " x ") %>% 
  filter(`Mating pair` != "A_IA_I x A_IA_I" &
           `Mating pair` != "A_OA_O x A_OA_O" &
           `Mating pair` != "X_IX_I x X_IY_O" &
           `Mating pair` != "X_OX_O x X_OY_O" &
           `Mating pair` != "Z_IW_O x Z_IZ_I" &
           `Mating pair` != "Z_OW_O x Z_OZ_O") %>% 
  mutate(Female = case_when(chromosome == "Autosome" | chromosome == "X" ~ 
                                  Female_homozygote_freq/
                                  (Female_homozygote_freq + Female_heterozygote_freq),
                                chromosome == "Z" ~ Female_hemizygote_freq),
         Male = case_when(chromosome == "Autosome" | chromosome == "Z" ~ 
                                  Male_homozygote_freq/
                                  (Male_homozygote_freq + Male_heterozygote_freq),
                                chromosome == "X" ~ Male_hemizygote_freq)) %>% 
  select(Generation, Dominance, `Mating pair`, expressed_sex, chromosome, Female,
         Male) %>% 
  pivot_longer(cols = Female:Male, names_to = "Sex", values_to = "prop_hom") %>% 
  mutate(prop_hom = if_else(prop_hom == "NaN", 0, prop_hom))

Autosomal loci

plotting_data %>% 
  filter(str_detect(`Mating pair`, "A")) %>% 
  ggplot(aes(x = Generation, 
                y = prop_hom,
                colour = Dominance, linetype = Sex)) +
  geom_line(linewidth = 1.75, alpha = 0.75) +
  scale_colour_manual(values = c("#38A6A5",
                                 "#73AF48",
                                 "#F89C74")) +
  labs(x = "Generation of inbreeding",
       y = "Prop. _I_ carriers that are homozygous",
       linetype =" Sex genotyped",
       title = "Autosomal loci") +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_continuous(expand = c(0,0),
                     limits = c(0, 20),
                     breaks = c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20)) +
  facet_grid(`Mating pair`~expressed_sex,
             labeller = label_glue(
      rows = '{`Mating pair`}',
      cols = 'Inbreeding controlled by {`expressed_sex`}')) +
  theme_bw() +
  theme(axis.title = element_markdown(),
        strip.background = element_rect(fill = "aliceblue"),
        plot.title = element_text(hjust = 0.5)) 

X-linked loci

plotting_data %>% 
  filter(str_detect(`Mating pair`, "X")) %>% 
  ggplot(aes(x = Generation, 
                y = prop_hom,
                colour = Dominance, linetype = Sex)) +
  geom_line(linewidth = 1.75, alpha = 0.75) +
  scale_colour_manual(values = c("#38A6A5",
                                 "#73AF48",
                                 "#F89C74")) +
  labs(x = "Generation of inbreeding",
       y = "Prop. _I_ carriers that are homozygous",
       linetype =" Sex genotyped",
       title = "X-linked loci") +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_continuous(expand = c(0,0),
                     limits = c(0, 20),
                     breaks = c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20)) +
  facet_grid(`Mating pair`~expressed_sex,
             labeller = label_glue(
      rows = '{`Mating pair`}',
      cols = 'Inbreeding controlled by {expressed_sex}')) +
  theme_bw() +
  theme(axis.title = element_markdown(),
        strip.background = element_rect(fill = "aliceblue"),
        plot.title = element_text(hjust = 0.5)) 

Z-linked loci

plotting_data %>% 
  filter(str_detect(`Mating pair`, "Z")) %>% 
  ggplot(aes(x = Generation, 
                y = prop_hom,
                colour = Dominance, linetype = Sex)) +
  geom_line(linewidth = 1.75, alpha = 0.75) +
  scale_colour_manual(values = c("#38A6A5",
                                 "#73AF48",
                                 "#F89C74")) +
  labs(x = "Generation of inbreeding",
       y = "Prop. _I_ carriers that are homozygous",
       linetype =" Sex genotyped",
       title = "Z-linked loci") +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_continuous(expand = c(0,0),
                     limits = c(0, 20),
                     breaks = c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20)) +
  facet_grid(`Mating pair`~expressed_sex,
             labeller = label_glue(
      rows = '{`Mating pair`}',
      cols = 'Inbreeding controlled by {expressed_sex}')) +
  theme_bw() +
  theme(axis.title = element_markdown(),
        strip.background = element_rect(fill = "aliceblue"),
        plot.title = element_text(hjust = 0.5)) 

The simulation

Build convenience functions

1. A sampling function that can handle vectors of length one

# so we can sample from vectors with length 1 without this being interpreted as an integer
sample_vec <- function(x, ...) x[sample(length(x), ...)] 

2. A function that builds the appropriate inheritance system.

make_mating_table <- function(gene_location){
  
  make_offspring <- function(X, Y, offspring_genotype, zygote_freq, sex){
    data.frame(Female_genotype = X,
           Male_genotype = Y,
           offspring_genotype,
           zygote_freq,
           sex)
  }
  
  # Specify the possible offspring genotypes for all the potential crosses; we use these for the offspring_genotype argument in the make_offspring function
  
  # offspring genotypes
  
  offspring_genotypes_1 <- c(2,2)
  offspring_genotypes_2 <- c(1, 1, 2, 2)
  offspring_genotypes_3 <- c(1, 1)
  offspring_genotypes_4 <- c(0, 0, 1, 1, 2, 2)
  offspring_genotypes_5 <- c(0, 0, 1, 1)
  offspring_genotypes_6 <- c(0, 0)
  
  offspring_genotypes_7 <- c(1, 2)
  offspring_genotypes_8 <- c(0, 1, 1, 2)
  offspring_genotypes_9 <- c(0, 1)
  offspring_genotypes_10 <- c(1, 0) # this is diff from above bc of the order with the sexes
  offspring_genotypes_11 <- c(2, 1)
  offspring_genotypes_12 <- c(2,1,1,0)
  
  # offspring sex
  
  offspring_sex_2 <- c(0, 1)
  offspring_sex_4 <- c(0, 1, 0, 1)
  offspring_sex_6 <- c(0, 1, 0, 1, 0, 1)
  
  # even frequency of two offspring genotypes
  
  freq_2 <- rep(0.5, 2)
  
  # even frequency between four offspring types
  
  freq_4 <- rep(0.25, 4)
  
  # when there are 6 offspring genotypes
  
  freq_6 <- c(0.125, 0.125,
              0.25, 0.25,
              0.125, 0.125)
  
  if(gene_location == "A"){
    books <- rbind(
        make_offspring(2, 2, offspring_genotypes_1, freq_2, offspring_sex_2),
        make_offspring(2, 1, offspring_genotypes_2, freq_4, offspring_sex_4),
        make_offspring(2, 0, offspring_genotypes_3, freq_2, offspring_sex_2),
        make_offspring(1, 2, offspring_genotypes_2, freq_4, offspring_sex_4),
        make_offspring(1, 1, offspring_genotypes_4, freq_6, offspring_sex_6),
        make_offspring(1, 0, offspring_genotypes_5, freq_4, offspring_sex_4),
        make_offspring(0, 2, offspring_genotypes_3, freq_2, offspring_sex_2),
        make_offspring(0, 1, offspring_genotypes_5, freq_4, offspring_sex_4),
        make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2)
    )
  }
  
  if(gene_location == "X"){
    books <- rbind(
        make_offspring(2, 1, offspring_genotypes_7, freq_2, offspring_sex_2),
        make_offspring(2, 0, offspring_genotypes_3, freq_2, offspring_sex_2),
        make_offspring(1, 1, offspring_genotypes_8, freq_4, offspring_sex_4),
        make_offspring(1, 0, offspring_genotypes_5, freq_4, offspring_sex_4),
        make_offspring(0, 1, offspring_genotypes_9, freq_2, offspring_sex_2),
        make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2)
    )
  }
  
  if(gene_location == "Y"){
    books <- rbind(
        make_offspring(0, 1, offspring_genotypes_10, freq_2, offspring_sex_2),
        make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2)
    )
  }
  
  if(gene_location == "Z"){
    books <- rbind(
        make_offspring(1, 2, offspring_genotypes_11, freq_2, offspring_sex_2),
        make_offspring(0, 2, offspring_genotypes_3, freq_2, offspring_sex_2),
        make_offspring(1, 1, offspring_genotypes_12, freq_4, offspring_sex_4),
        make_offspring(0, 1, offspring_genotypes_5, freq_4, offspring_sex_4),
        make_offspring(1, 0, offspring_genotypes_10, freq_2, offspring_sex_2),
        make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2)
    )
  }
  
  if(gene_location == "W"){
    books <- rbind(
        make_offspring(1, 0, offspring_genotypes_9, freq_2, offspring_sex_2),
        make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2)
    )
  }
  
  if(gene_location == "C"){
    books <- rbind(
      make_offspring(1, 0, offspring_genotypes_3, freq_2, offspring_sex_2),
      make_offspring(1, 1, offspring_genotypes_3, freq_2, offspring_sex_2),
      make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2),
      make_offspring(0, 1, offspring_genotypes_6, freq_2, offspring_sex_2)
    )
  }
  
  if(gene_location == "P"){
    books <- rbind(
      make_offspring(1, 0, offspring_genotypes_6, freq_2, offspring_sex_2),
      make_offspring(1, 1, offspring_genotypes_3, freq_2, offspring_sex_2),
      make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2),
      make_offspring(0, 1, offspring_genotypes_3, freq_2, offspring_sex_2)
    )
  }
    return(books)  
}


offspring_genotypes_autosome <- make_mating_table("A")
offspring_genotypes_X <- make_mating_table("X")
offspring_genotypes_Y <- make_mating_table("Y")
offspring_genotypes_Z <- make_mating_table("Z")
offspring_genotypes_W <- make_mating_table("W")
offspring_genotypes_C <- make_mating_table("C")
offspring_genotypes_P <- make_mating_table("P")

3. A function that takes two parental genotypes and produces offspring

sample_mating_table <- function(inheritance_scheme, 
                                f,
                                mother){
  
  # cut to possible genotypes
  possibilities <- 
    inheritance_scheme[inheritance_scheme$Female_genotype == mother[4] &
                         inheritance_scheme$Male_genotype == mother[9], c(3,5)]
  # get prob of producing each genotype
  probs <- 
    inheritance_scheme[inheritance_scheme$Female_genotype == mother[4] &
                         inheritance_scheme$Male_genotype == mother[9], 4]
  # sample
  possibilities[sample(size = f,
                       x = nrow(possibilities), 
                       prob = probs,
                       replace = TRUE), ]
}

Load the parameter space

resolution <- 25
starting_pop_size_autosomes <- 2000 # both sexes harbour two copies of each autosomal chromosome = 1000 autosomal haplotypes

parameters <-
  expand_grid(
    chromosome = c("A", "X", "Y", "Z", "W", "C", "P"),
    v = c(8, 80),
    D = seq(0, -0.99, length = resolution), # inbreeding depression
    refractory_period_prop_cohort_alive = seq(0.01, 1, length = resolution)
  ) %>% 
  full_join(tibble(chromosome = c("A", "A", "A", "A", "A", "A",
                                  "X", "X", "X", "X",
                                  "Y",
                                  "Z", "Z", "Z", "Z",
                                  "W",
                                  "C", "C",
                                  "P", "P"),
                   sex_expressed = c(0, 0, 0, 1, 1, 1,
                                     0, 1, 1, 1,
                                     0,
                                     0, 0, 0, 1,
                                     1,
                                     0, 1,
                                     0, 1),
                   dominance = c(0, 0.5, 1, 0, 0.5, 1,
                                 1, 0, 0.5, 1, 
                                 1, 
                                 0, 0.5, 1, 1,
                                 1,
                                 1, 1,
                                 1, 1)) %>% 
              mutate(Starting_pop_size = case_when(chromosome == "A" ~ starting_pop_size_autosomes,
                                                   chromosome == "X" | chromosome == "Z" ~ 
                                                     starting_pop_size_autosomes / 0.75,
                                                   chromosome == "Y" | chromosome == "W" ~ starting_pop_size_autosomes*4,
                                                   chromosome == "C" | chromosome == "P" ~ starting_pop_size_autosomes*2)),
            relationship = "many-to-many", by = "chromosome") %>% 
  mutate(baseline_mean_lifespan = 1,
         v = v / (Starting_pop_size / 2),
         f = 5, 
         refractory_period = -log(refractory_period_prop_cohort_alive),
         mutation_time = 5, # this is when the mutation can be introduced from
         time_end = 1000, # with avg lifespan = 1, this is ~ roughly 1000 gens
         parameter_space_ID = row_number(),
         mutation_events = 5)

parameters_autosome <- parameters %>% filter(chromosome == "A") %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs
parameters_X <- parameters %>% filter(chromosome == "X") %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs
parameters_Y <- parameters %>% filter(chromosome == "Y") %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs
parameters_Z <- parameters %>% filter(chromosome == "Z") %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs
parameters_W <- parameters %>% filter(chromosome == "W") %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs
parameters_C <- parameters %>% filter(chromosome == "C") %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs
parameters_P <- parameters %>% filter(chromosome == "P") %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs

if(!file.exists("parameters/parameters_autosome.txt")){
  parameters_autosome %>% write.table("parameters/parameters_autosome.txt")
  parameters_X %>% write.table("parameters/parameters_X.txt")
  parameters_Y %>% write.table("parameters/parameters_Y.txt")
  parameters_Z %>% write.table("parameters/parameters_Z.txt")
  parameters_W %>% write.table("parameters/parameters_W.txt")
  parameters_C %>% write.table("parameters/parameters_C.txt")
  parameters_P %>% write.table("parameters/parameters_P.txt")
}

The simulation function

continuous_time_simulation <- function(row,
                                       parameters,
                                       inheritance_scheme){
  
  #print(paste("Doing row", row)) # this shows which row in the parameter space is being modelled
  
  Starting_pop_size <- round(parameters$Starting_pop_size[row], 0)
  f <- parameters$f[row] # fecundity constant
  mutation_time <- parameters$mutation_time[row] # introduce an I allele after family structure is established
  baseline_mean_lifespan <- parameters$baseline_mean_lifespan[row] # constant at 1
  time_end <- parameters$time_end[row] # a cut-off point for each run 
  sex_expressed <- parameters$sex_expressed[row]
  chromosome <- parameters$chromosome[row]
  v <- parameters$v[row]
  refractory_period <- parameters$refractory_period[row]
  D <- parameters$D[row]
  dominance <- parameters$dominance[row]
  parameter_space_ID <- parameters$parameter_space_ID[row]
  mutation_events <- parameters$mutation_events[row]
  
  # Set the number of breeding sites
  
  breeding_sites <- round(0.2*Starting_pop_size, 0)
  
  # what inheritance system does this run follow
  offspring_genotypes <- inheritance_scheme
  
  # Set the maximum number of I alleles that can be found in each sex
  if(chromosome == "A"){
    female_max_I <- 2
    male_max_I <- 2
  }
  
  if(chromosome == "X"){
    female_max_I <- 2
    male_max_I <- 1
  }
  
  if(chromosome == "Z"){
    female_max_I <- 1
    male_max_I <- 2
  }
  
  if(chromosome == "Y"){
    female_max_I <- 0
    male_max_I <- 1
  }
  
  if(chromosome == "W"){
    female_max_I <- 1
    male_max_I <- 0
  }
  
  if(chromosome == "C" | chromosome == "P"){
    female_max_I <- 1
    male_max_I <- 1
  }
  
  # make matrix to hold results; updated as sim progresses
  # col1 = time, col2 = prop I, col3 = pop size, col4 = prop virgin female deaths
  results_matrix <- matrix(nrow = time_end*4+2, ncol = 4) # record each time point
  
  # make matrix to hold population; updated as sim progresses
  
  # col1 = ID 
  # col2 = Family ID
  # col3 = Sex: females = 1 and males = 0
  # col4 = Genotype: 0, 1 and 2 = copies of inbreeding allele
  # col5 = mortality rate
  # col6 = encountered relative: NA = NO, 1 = YES
  # col7 = mating state: -Inf not in pop, NA = unmated, real = out, Inf = mated female
  # col8 = inbred mating: NA = NO, 1 = YES (only matters for females)
  # col9 = mated_genotype: NA = unmated, otherwise 0,1,2 (see mating table)
  # col10 = breeding site:  NA = NO, 1 = YES 
  # col11 = no. matings (only matters for males)
  # col12 = offspring produced: NA = NO, 1 = YES
  
  pop_matrix <- matrix(nrow = Starting_pop_size*2, # pop expands with initial repro pulse
                       ncol = 12)
  # ID & Family ID
  pop_matrix[1:Starting_pop_size, 1:2] <- 1:Starting_pop_size
  # assign sex
  pop_matrix[1:Starting_pop_size, 3] <- rbinom(n = Starting_pop_size, 1, prob = 0.5)
  # female_starting_genotype
  pop_matrix[pop_matrix[,3] < 1 & !is.na(pop_matrix[,3]), 4] <- 0
  # male_starting_genotype
  pop_matrix[pop_matrix[,3] > 0 & !is.na(pop_matrix[,3]), 4] <- 0
  # assign mortality rates
  pop_matrix[1:Starting_pop_size, 5] <- 1/baseline_mean_lifespan
  # set the unused rows to state -Inf 
  pop_matrix[(Starting_pop_size + 1):nrow(pop_matrix), 7] <- -Inf
  # mate count
  pop_matrix[1:Starting_pop_size, 11] <- 0
  # offspring production status
  pop_matrix[1:Starting_pop_size, 12] <- NA
  # populate breeding sites
  # the starting no. of females generally exceeds the number of breeding sites, which is starting_pop_size/f. The code below selects the initial breeding site holders
  
  if(nrow(pop_matrix[pop_matrix[,3] > 0 & !is.na(pop_matrix[,3]),]) > breeding_sites){
    
    initial_breeders <- head(pop_matrix[pop_matrix[,3] > 0 & !is.na(pop_matrix[,3]),1], breeding_sites)
    
    pop_matrix[initial_breeders,10] <- 1 # take advantage of ID = row number for initial pop
    
  } else{pop_matrix[pop_matrix[,3] > 0 & !is.na(pop_matrix[,3]), 10] <- 1}
  
  # Initialise counter for the results table
  
  next_update <- 0 # keep track of when to update the results
  next_row <- 0 # keep track of which row to update
  
  # Initialise the Individual_ID and Family_ID counters
  
  Individual_ID_counter <- Starting_pop_size
  
  Family_ID_counter <- Starting_pop_size # each individual descends from a distinct family at onset
  
  # Initialise the timer t
  
  t <- 0
  
  # Set initial pop size and freq of I allele for results table
  
  Prop_I <- 0 
  pop_size <- Starting_pop_size
  total_female_deaths <- 0
  mated_female_deaths <- 0
  
  # Start population without the I allele to generate family structure
  # Flips to 1 at mutant intro time point 
  
  mutant_introduced <- 0
  
  keep_going <- TRUE # if the inbreeding allele fixes or goes extinct, this will change to false and the while loop will quit early
  
  # With the initial population ready to go, start the timer and let the simulation run.
  
  while(t <= time_end & keep_going){
    
    print(paste0("Population size = ", pop_size, 
                 ", breeders = ", sum(pop_matrix[,10], na.rm = T), 
                 ", time = ", round(t, 3), ", Prop I =", Prop_I, ", mutation events =", mutant_introduced))
    
    # find next event 
    
    # next death: this is the sum of the mortality rates for all individuals in the population
    
    next_death <- t + rexp(n = 1, rate = sum(pop_matrix[, 5], na.rm = T))
    
    # next receptive mating encounter
    
    # find no. of females in mating pool & separate by encounter experience
    
    receptive_females_first_encounter <- 
      pop_matrix[pop_matrix[,3] > 0 &
                   is.na(pop_matrix[,6]) &
                   is.na(pop_matrix[,7]),, drop = FALSE]
    
    receptive_females_second_encounter <- 
      pop_matrix[pop_matrix[,3] > 0 &
                   !is.na(pop_matrix[,6]) &
                   is.na(pop_matrix[,7]),, drop = FALSE]
    
    # find no. of males in mating pool
    receptive_males <- pop_matrix[pop_matrix[,3] < 1 & is.na(pop_matrix[,7]),, drop = FALSE]
    
    # Find the time the next encounter occurs: plug the sum of the rates into the exponential function. 
    # The population level encounter rate is the product of the rate at which a single male finds a single female, the number of receptive females in the population, and the number of receptive males in the population
    
    next_first_encounter <- t + 
      rexp(n = 1, rate = v*nrow(receptive_females_first_encounter)*nrow(receptive_males))
    
    next_secondary_encounter <- t + 
      rexp(n = 1, rate = v*nrow(receptive_females_second_encounter)*nrow(receptive_males))
    
    # time in - Inf, Inf and Na are possible options that the code can handle 
    next_time_in <- min(pop_matrix[is.finite(pop_matrix[,7]),7])
    
    # find which event happens next and update t
    t <- pmin(next_death,
              next_time_in, 
              next_first_encounter,
              next_secondary_encounter,
              next_update, # update the population
              na.rm = TRUE) # ... if a rate is 0, NaN produced.
    
    
    if(t == next_update & !is.na(next_update)){# record time, I prop and pop size
      results_matrix[next_row+1,1] <- t
      results_matrix[next_row+1,2] <- round(Prop_I, 4)
      results_matrix[next_row+1,3] <- pop_size # popsize
      results_matrix[next_row+1,4] <- round(mated_female_deaths / total_female_deaths, 3)
      next_update <- next_update + 0.25
      next_row <- next_row + 1
      total_female_deaths <- 0 # reset the count
      mated_female_deaths <- 0 # reset the count
    }
    
    
    if(t == next_death){# remove an individual from the pop
      who_died <- 
        sample_vec(size = 1, # choose one
                   x = pop_matrix[!is.na(pop_matrix[,1]),1], # subset to current pop
                   prob = pop_matrix[!is.na(pop_matrix[,5]),5]) # weight by mortality rate
      # add a death if it was a female
      if(nrow(pop_matrix[pop_matrix[,1] == who_died &
                   !is.na(pop_matrix[,1]) &
                   pop_matrix[,3] > 0,, drop = FALSE]) > 0){total_female_deaths <- total_female_deaths + 1}
      
      # add virgin female deaths
      if(nrow(pop_matrix[pop_matrix[,1] == who_died &
                   !is.na(pop_matrix[,1]) &
                   pop_matrix[,3] > 0 &
                   is.infinite(pop_matrix[,7]),, drop = FALSE]) > 0){mated_female_deaths <- mated_female_deaths + 1}
      
      # remove individual from pop matrix
      pop_matrix[pop_matrix[,1] == who_died, 7] <- -Inf # NA means time-in here, so special edit required
      pop_matrix[pop_matrix[,1] == who_died, c(1:6, 8:12)] <- NA 
      
      # re-order to make steps like adding offspring easier later on
      pop_matrix <- pop_matrix[order(pop_matrix[,1]),]
      
    }
    
    # check if there are free breeding sites and whether females are available to fill them 
    
    current_breeders <- sum(pop_matrix[, 10], na.rm = T)
    
    # get list of IDs for floating females
    floating_females <- pop_matrix[!is.na(pop_matrix[,1]) & # alive
                                     pop_matrix[,3] > 0 & # female
                                     is.na(pop_matrix[,10]), # non-breeding
                                   1] # return the IDs only
    
    # If so, recruit a new breeder
    # All prospective females have equal probability
    
    if(current_breeders < breeding_sites & length(floating_females) > 0){
      
      # assign the new breeder
      
      new_breeder <- 
        sample_vec(size = 1, # choose one
                   x = floating_females) # subset to floaters
      
      pop_matrix[pop_matrix[,1] == new_breeder, 10] <- 1
    }
    
    if(t == next_time_in & !is.na(next_time_in)){ # a male re-enters the mating pool
      pop_matrix[pop_matrix[,7] == next_time_in, 7] <- NA # change to receptive
    }
    
    #### mating
    
    if(t == next_first_encounter &
       !is.na(next_first_encounter)){# does first encounter lead to (inbred) mating?
      
      # Determine whether a heterozygote inbreeds on this occasion. 
      # Depends on genotype if this matters
      heterozygote_inbreeds <- rbinom(1, 1, prob = dominance)
      
      # which female
      female_ID <- sample_vec(receptive_females_first_encounter[,1], 1)
      # get meta-data
      female <- subset(pop_matrix, pop_matrix[,1] == female_ID)
      # how many inbreeding alleles does she carry?
      alleles_female <- female[,4]
      
      mates <- NULL # reset this every time as a safeguard - MAYBE REMOVE?
      
      # find brothers that are in the mating pool
      brothers <-
        pop_matrix[pop_matrix[,2] == female[, 2] & # find family members
                     pop_matrix[,3] < 1 & # that are male
                     is.na(pop_matrix[,7]) & # and in the mating pool
                     !is.na(pop_matrix[,1]), # remove NAs
                   1] 
      # find the specific brother - if there aren't any, inbreeding does not happen
      if(length(brothers) > 0){# choose brother randomly
        chosen_brother <-
          subset(pop_matrix, 
                 pop_matrix[,1] == sample_vec(size = 1, x = brothers))
        # how many inbreeding alleles does he carry?
        alleles_brother <- chosen_brother[,4]
        brother_ID <- chosen_brother[,1]
      }else{alleles_brother <- 0} # we need this for the next if statement
      
      # now determine whether inbreeding occurs:
      # which individual expresses the allele
      # does that individual have the allele
      # is it expressed (depends on genomic region, no. copies and dominance)
      
      if(# female expression determines outcome
        # dominance doesn't matter
        length(brothers) > 0 & sex_expressed > 0 & female_max_I == alleles_female |
        # dominance matters
        length(brothers) > 0 & sex_expressed > 0 & 
        0 < alleles_female & alleles_female < female_max_I & heterozygote_inbreeds > 0 |
        # male expression determines outcome
        # dominance doesn't matter
        length(brothers) > 0 & sex_expressed < 1 & male_max_I == alleles_brother |
        # dominance matters
        length(brothers) > 0 & sex_expressed < 1 & 
        0 < alleles_brother & alleles_brother < male_max_I & heterozygote_inbreeds > 0){
        
        # do inbreeding
        # update the pop matrix
        # female
        pop_matrix[pop_matrix[,1] == female_ID, 6] <- 1 # relative has been encountered
        pop_matrix[pop_matrix[,1] == female_ID, 7] <- Inf # female leaves mating pool
        pop_matrix[pop_matrix[,1] == female_ID, 8] <- 1 # inbreeding occurs
        pop_matrix[pop_matrix[,1] == female_ID, 9] <- alleles_brother # mates genotype
        
        # male
        pop_matrix[pop_matrix[,1] == brother_ID, 7] <- t + refractory_period # male leaves mating pool
        pop_matrix[pop_matrix[,1] == brother_ID, 8] <- 1 # inbreeding occurs
        pop_matrix[pop_matrix[,1] == brother_ID & !is.na(pop_matrix[,1]), 11] <-
          pop_matrix[pop_matrix[,1] == brother_ID & !is.na(pop_matrix[,1]), 11] + 1
      } else{
        # inbreeding is avoided
        # females that had no receptive brother to encounter are recorded as having had their chance for inbreeding early in life. When the male refractory period != 0, this is possible but unlikely (because all siblings are produced at the same time). Most commonly, this will occur when a female produces an all-female brood (0.03125 probability when f=5)
        
        pop_matrix[pop_matrix[,1] == female_ID, 6] <- 1 # relative has been encountered
      }
    }
    
    if(t == next_secondary_encounter &
       !is.na(next_secondary_encounter)){ 
      # If the individual has already encountered a sibling, don't swap and let encounter proceed. 
      
      # which female
      female_ID <- sample_vec(receptive_females_second_encounter[,1], 1)
      # get meta-data
      female <- subset(pop_matrix, pop_matrix[,1] == female_ID)
      # how many inbreeding alleles does she carry?
      alleles_female <- female[,4]  
      
      # which male
      male_ID <- sample_vec(receptive_males[,1], 1)
      # get meta-data
      male <- subset(pop_matrix, pop_matrix[,1] == male_ID)
      # how many inbreeding alleles does he carry?
      alleles_male <- male[,4] 
      
      # If the pair happen to be siblings, check if they inbreed  
      
      # Determine whether a heterozygote inbreeds on this occasion. 
      # Depends on genotype if this matters
      heterozygote_inbreeds <- rbinom(1, 1, prob = dominance)
      
      if(
        # female expression determines outcome
        # dominance doesn't matter
        female[,2] == male[,2] & sex_expressed > 0 & female_max_I == alleles_female |
        # dominance matters
        female[,2] == male[,2] & sex_expressed > 0 & 
        0 < alleles_female & alleles_female < female_max_I & heterozygote_inbreeds > 0 |
        # male expression determines outcome
        # dominance doesn't matter
        female[,2] == male[,2] & sex_expressed < 1 & male_max_I == alleles_male |
        # dominance matters
        female[,2] == male[,2] & sex_expressed < 1 & 
        0 < alleles_male & alleles_male < male_max_I & heterozygote_inbreeds > 0){
        
        # do inbreeding
        # update the pop matrix
        # female
        pop_matrix[pop_matrix[,1] == female_ID, 7] <- Inf # female leaves mating pool
        pop_matrix[pop_matrix[,1] == female_ID, 8] <- 1 # inbreeding occurs
        pop_matrix[pop_matrix[,1] == female_ID, 9] <- alleles_male # mates genotype
        
        # male
        pop_matrix[pop_matrix[,1] == male_ID, 7] <- t + refractory_period # male leaves mating pool
        pop_matrix[pop_matrix[,1] == male_ID & !is.na(pop_matrix[,1]), 11] <-
          pop_matrix[pop_matrix[,1] == male_ID & !is.na(pop_matrix[,1]), 11] + 1
      } else{
        # do outbreeding
        # update the pop matrix
        # female
        pop_matrix[pop_matrix[,1] == female_ID, 7] <- Inf # female leaves mating pool
        pop_matrix[pop_matrix[,1] == female_ID, 9] <- alleles_male # mates genotype
        
        # male
        pop_matrix[pop_matrix[,1] == male_ID, 7] <- t + refractory_period # male leaves mating pool
        pop_matrix[pop_matrix[,1] == male_ID & !is.na(pop_matrix[,1]), 11] <-
          pop_matrix[pop_matrix[,1] == male_ID & !is.na(pop_matrix[,1]), 11] + 1
      }
    }
    
    # Consequences of death and mating: reproduction
    
    # check if a female can now produce offspring, either because they're previously mated and have secured a breeding site or because they already hold a breeding site and have now mated
    # make sure that previous breeders are excluded
    
    new_mated_breeder <- pop_matrix[is.infinite(pop_matrix[,7]) & # mated
                                      !is.na(pop_matrix[,10]) & # holds breeding site
                                      is.na(pop_matrix[,12]),, drop = FALSE] # hasn't reproduced
    
    if(nrow(new_mated_breeder) > 0){
      # add offspring to the population
      # each mated female that holds a breeding site produces f offspring
      
      # first check whether the mutant I allele should be added
      if(mutant_introduced < mutation_events & t > mutation_time){
        which_sex <- rbinom(1, 1, prob = 0.5)
        
        if(chromosome == "A" & which_sex == 1 |
           chromosome == "X" & which_sex == 1 |
           chromosome == "Z" & which_sex == 1){
          new_mated_breeder[4] <- 1
        }
        
        if(chromosome == "A" & which_sex == 0 |
           chromosome == "X" & which_sex == 0 |
           chromosome == "Z" & which_sex == 0){
          new_mated_breeder[9] <- 1
        }
        
        if(chromosome == "W"|
           chromosome == "C"){
          new_mated_breeder[4] <- 1
        }
        
        if(chromosome == "Y" |
           chromosome == "P"){
          new_mated_breeder[9] <- 1
        }
        
        mutant_introduced <- mutant_introduced + 1
      }
      
      next_row_to_fill <- length(pop_matrix[!is.na(pop_matrix[,1]),1]) + 1
      last_row_to_fill <- next_row_to_fill + f - 1
      next_ID <- Individual_ID_counter + 1
      last_ID <- Individual_ID_counter + f
      Family_ID_counter <- Family_ID_counter + 1
      
      # assign IDs
      pop_matrix[next_row_to_fill:last_row_to_fill, 1] <- next_ID:last_ID
      # assign all offspring to a single family
      pop_matrix[next_row_to_fill:last_row_to_fill, 2] <- Family_ID_counter
      # assign sex and genotype using our mating table sampling function
      offspring_genos <- 
        sample_mating_table(inheritance_scheme,
                            f, 
                            mother = new_mated_breeder)
      pop_matrix[next_row_to_fill:last_row_to_fill, 3] <- offspring_genos[,2]
      pop_matrix[next_row_to_fill:last_row_to_fill, 4] <- offspring_genos[,1]
      # assign mortality rates
      if(is.na(new_mated_breeder[8])){
        pop_matrix[next_row_to_fill:last_row_to_fill, 5] <- 1/baseline_mean_lifespan
      } else{ # apply effect of inbreeding depression
        pop_matrix[next_row_to_fill:last_row_to_fill, 5] <- 1/(baseline_mean_lifespan + D)
      }
      # fill in the mating and breeding site details - everyone starts as a floating virgin
      pop_matrix[next_row_to_fill:last_row_to_fill, 6:10] <- NA
      # mate count
      pop_matrix[next_row_to_fill:last_row_to_fill, 11] <- 0
      
      
      # update the mothers offspring production status
      
      pop_matrix[pop_matrix[,1] == new_mated_breeder[1], 12] <- 1
      
      # update the individual ID counter (redundant but more readable to do this here)
      Individual_ID_counter <- last_ID
      
    }      
    
    # Calculate the frequency of the I allele, quit early if I fixes or goes extinct
    
    pop_size <- nrow(pop_matrix[!is.na(pop_matrix[,1]),, drop = FALSE]) # use this to update the results
    n_females <- nrow(pop_matrix[!is.na(pop_matrix[,1]) &
                                   pop_matrix[,3] > 0,, drop = FALSE])
    n_males <- pop_size - n_females
    
    # calc allele freq if autosomal locus   
    if(chromosome == "A"){
      Prop_I <-
        sum(pop_matrix[,4], na.rm = T)/(pop_size*2) # x2 because diploid
    }
    
    # calc allele freq if W locus   
    if(chromosome == "W"){
      Prop_I <-
        sum(pop_matrix[,4], na.rm = T)/n_females
    }
    
    # calc allele freq if Y locus   
    if(chromosome == "Y"){
      Prop_I <-
        sum(pop_matrix[,4], na.rm = T)/n_males 
    }
    
    # calc allele freq if X locus   
    if(chromosome == "X"){
      Prop_I <-
        sum(pop_matrix[,4], na.rm = T)/(n_females*2 + n_males)
    }
    
    # calc allele freq if Z locus   
    if(chromosome == "Z"){
      Prop_I <-
        sum(pop_matrix[,4], na.rm = T)/(n_females + n_males*2) 
    }
    
    # calc allele freq if C locus   
    if(chromosome == "C" |
       chromosome == "P"){
      Prop_I <-
        sum(pop_matrix[,4], na.rm = T)/pop_size 
    }
    
    # quit condition
    if(mutant_introduced > 0 & Prop_I > 0.9 |
       mutant_introduced > 0 & Prop_I == 0 | 
       pop_size < 2){keep_going <- FALSE}
    
  }
  
  results_matrix[next_row+1,1] <- t
  results_matrix[next_row+1,2] <- round(Prop_I, 4)
  results_matrix[next_row+1,3] <- pop_size
  results_matrix[next_row+1,4] <- round(mated_female_deaths / total_female_deaths, 3)
  results_matrix <- results_matrix[-(next_row+2:nrow(results_matrix)),]
  # save results as a csv.  
  
  results_matrix
  
  write.csv(results_matrix,
            paste("results/rowID_", 
                  parameter_space_ID, 
                  chromosome, ".csv", 
                  sep = ""))
  #write.csv(results_matrix,
    #        paste("sim_results/rowID_", 
      #            parameter_space_ID, 
        #          chromosome, ".csv", 
          #        sep = ""))
}

Run the simulation

In practice, I ran the simulations on JGU’s Mogon computing cluster. See the HPC_inbreeding_script.R and the batch script run_inbreeding_sim. To run the simulation for a single parameter space, you could run continuous_time_simulation(1, parameters_P, offspring_genotypes_P)

Load the results

# build a function to load the individual runs and join them into a single tibble

files <-
    list.files(path = "sim_results") %>% 
    str_subset("P") # change this to load the desired files

results_reader <- function(x){
  read_csv(paste0("sim_results/", x)) %>% 
    mutate(parameter_space_ID = x)
}

if(!file.exists("results/autosome_results.csv")){
  results <- 
    map_dfr(files, results_reader) %>% 
    rename(time = V1,
           I_prop = V2,
           pop_size = V3,
           prop_mated = V4) %>% 
    select(-`...1`) %>% 
    mutate(parameter_space_ID = str_remove(parameter_space_ID, "rowID_"),
           parameter_space_ID = as.integer(str_remove(parameter_space_ID, "A.csv"))) %>% 
    left_join(parameters_autosome, by = "parameter_space_ID")
  
  write_csv(results, "results/autosome_results.csv")
}else{
  autosome_results <- read_delim("results/autosome_results.csv")
}

if(!file.exists("results/X_results.csv")){
  results <- 
    map_dfr(files, results_reader) %>% 
    rename(time = V1,
           I_prop = V2,
           pop_size = V3,
           prop_mated = V4) %>% 
    select(-`...1`) %>% 
    mutate(parameter_space_ID = str_remove(parameter_space_ID, "rowID_"),
           parameter_space_ID = as.integer(str_remove(parameter_space_ID, "X.csv"))) %>% 
    left_join(parameters_X, by = "parameter_space_ID")
  
  write_csv(results, "results/X_results.csv")
}else{
  X_results <- read_delim("results/X_results.csv")
}

if(!file.exists("results/Y_results.csv")){
  results <- 
    map_dfr(files, results_reader) %>% 
    rename(time = V1,
           I_prop = V2,
           pop_size = V3,
           prop_mated = V4) %>% 
    select(-`...1`) %>% 
    mutate(parameter_space_ID = str_remove(parameter_space_ID, "rowID_"),
           parameter_space_ID = as.integer(str_remove(parameter_space_ID, "Y.csv"))) %>% 
    left_join(parameters_Y, by = "parameter_space_ID")
  
  write_csv(results, "results/Y_results.csv")
}else{
  Y_results <- read_delim("results/Y_results.csv")
}

if(!file.exists("results/Z_results.csv")){
  results <- 
    map_dfr(files, results_reader) %>% 
    rename(time = V1,
           I_prop = V2,
           pop_size = V3,
           prop_mated = V4) %>% 
    select(-`...1`) %>% 
    mutate(parameter_space_ID = str_remove(parameter_space_ID, "rowID_"),
           parameter_space_ID = as.integer(str_remove(parameter_space_ID, "Z.csv"))) %>% 
    left_join(parameters_Z, by = "parameter_space_ID")
  
  write_csv(results, "results/Z_results.csv")
}else{
  Z_results <- read_delim("results/Z_results.csv")
}

if(!file.exists("results/W_results.csv")){
  results <- 
    map_dfr(files, results_reader) %>% 
    rename(time = V1,
           I_prop = V2,
           pop_size = V3,
           prop_mated = V4) %>% 
    select(-`...1`) %>% 
    mutate(parameter_space_ID = str_remove(parameter_space_ID, "rowID_"),
           parameter_space_ID = as.integer(str_remove(parameter_space_ID, "W.csv"))) %>% 
    left_join(parameters_W, by = "parameter_space_ID")
  
  write_csv(results, "results/W_results.csv")
}else{
  W_results <- read_delim("results/W_results.csv")
}

if(!file.exists("results/C_results.csv")){
  results <- 
    map_dfr(files, results_reader) %>% 
    rename(time = V1,
           I_prop = V2,
           pop_size = V3,
           prop_mated = V4) %>% 
    select(-`...1`) %>% 
    mutate(parameter_space_ID = str_remove(parameter_space_ID, "rowID_"),
           parameter_space_ID = as.integer(str_remove(parameter_space_ID, "C.csv"))) %>% 
    left_join(parameters_C, by = "parameter_space_ID")
  
  write_csv(results, "results/C_results.csv")
}else{
  C_results <- read_delim("results/C_results.csv")
}

if(!file.exists("results/P_results.csv")){
  results <- 
    map_dfr(files, results_reader) %>% 
    rename(time = V1,
           I_prop = V2,
           pop_size = V3,
           prop_mated = V4) %>% 
    select(-`...1`) %>% 
    mutate(parameter_space_ID = str_remove(parameter_space_ID, "rowID_"),
           parameter_space_ID = as.integer(str_remove(parameter_space_ID, "P.csv"))) %>% 
    left_join(parameters_P, by = "parameter_space_ID")
  
  write_csv(results, "results/P_results.csv")
}else{
  P_results <- read_delim("results/P_results.csv")
}

Plotting

temp <- pnw_palette("Shuksan2",100)

Population dynamics

Autosomes

a1 <-
  autosome_results %>%
  filter(chromosome == "A" & !(time %% 1)) %>%   
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.008 ~ "low",
                           .default = "high")) %>% 
  ggplot(aes(x = time, y = pop_size, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
        panel.spacing.x = unit(6, "mm"))

a2 <-
autosome_results %>%
  filter(chromosome == "A" & !(time %% 1)) %>% 
   mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.008 ~ "low",
                           .default = "high")) %>%  
  ggplot(aes(x = time, y = prop_mated, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Prop. mated",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
        panel.spacing.x = unit(6, "mm"))

a1 / a2 + plot_layout(guides = "collect") & theme(legend.position = 'bottom')

X chromosome

X1 <-
  X_results %>%
  filter(chromosome == "X" & !(time %% 1)) %>%   
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.006 ~ "low",
                           .default = "high")) %>% 
  ggplot(aes(x = time, y = pop_size, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 201)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
         panel.spacing.x = unit(6, "mm"))

X2 <-
X_results %>%
  filter(chromosome == "X" & !(time %% 1)) %>% 
   mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.006 ~ "low",
                           .default = "high")) %>%  
  ggplot(aes(x = time, y = prop_mated, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 201)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Prop. mated",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
         panel.spacing.x = unit(6, "mm"))

X1 + X2 + plot_layout(guides = "collect") & theme(legend.position = 'bottom')

Z chromosome

Z1 <-
  Z_results %>%
  filter(chromosome == "Z" & !(time %% 1)) %>%   
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.006 ~ "low",
                           .default = "high")) %>% 
  ggplot(aes(x = time, y = pop_size, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
         panel.spacing.x = unit(6, "mm"))

Z2 <-
Z_results %>%
  filter(chromosome == "Z" & !(time %% 1)) %>% 
   mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.006 ~ "low",
                           .default = "high")) %>%  
  ggplot(aes(x = time, y = prop_mated, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Prop. mated",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
         panel.spacing.x = unit(6, "mm"))

Z1 + Z2 + plot_layout(guides = "collect") & theme(legend.position = 'bottom')

W chromosome

W1 <-
  W_results %>%
  filter(chromosome == "W" & !(time %% 1)) %>%   
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.002 ~ "low",
                           .default = "high")) %>% 
  ggplot(aes(x = time, y = pop_size, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
         panel.spacing.x = unit(6, "mm"))

W2 <-
W_results %>%
  filter(chromosome == "W" & !(time %% 1)) %>% 
   mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.002 ~ "low",
                           .default = "high")) %>%  
  ggplot(aes(x = time, y = prop_mated, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Prop. mated",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
         panel.spacing.x = unit(6, "mm"))

W1 + W2 + plot_layout(guides = "collect") & theme(legend.position = 'bottom')

Y chromosome

Y1 <-
  Y_results %>%
  filter(chromosome == "Y" & !(time %% 1)) %>%   
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.002 ~ "low",
                           .default = "high")) %>% 
  ggplot(aes(x = time, y = pop_size, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
         panel.spacing.x = unit(6, "mm"))

Y2 <-
Y_results %>%
  filter(chromosome == "Y" & !(time %% 1)) %>% 
   mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.002 ~ "low",
                           .default = "high")) %>%  
  ggplot(aes(x = time, y = prop_mated, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Prop. mated",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
         panel.spacing.x = unit(6, "mm"))

Y1 + Y2 + plot_layout(guides = "collect") & theme(legend.position = 'bottom')

Maternally inherited cytoplasmic chromosome

C1 <-
  C_results %>%
  filter(chromosome == "C" & !(time %% 1)) %>%   
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.004 ~ "low",
                           .default = "high")) %>% 
  ggplot(aes(x = time, y = pop_size, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
         panel.spacing.x = unit(6, "mm"))

C2 <-
C_results %>%
  filter(chromosome == "C" & !(time %% 1)) %>% 
   mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.004 ~ "low",
                           .default = "high")) %>%  
  ggplot(aes(x = time, y = prop_mated, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Prop. mated",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
         panel.spacing.x = unit(6, "mm"))

C1 + C2 + plot_layout(guides = "collect") & theme(legend.position = 'bottom')

Paternally inherited cytoplasmic chromosome

P1 <-
  P_results %>%
  filter(chromosome == "P" & !(time %% 1)) %>%   
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.004 ~ "low",
                           .default = "high")) %>% 
  ggplot(aes(x = time, y = pop_size, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
         panel.spacing.x = unit(6, "mm"))

P2 <-
P_results %>%
  filter(chromosome == "P" & !(time %% 1)) %>% 
   mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.004 ~ "low",
                           .default = "high")) %>%  
  ggplot(aes(x = time, y = prop_mated, colour = 1-refractory_period_prop_cohort_alive)) + 
  geom_line(aes(group = parameter_space_ID), alpha = 0.6) +
  geom_vline(xintercept = 10, linetype = 2) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Prop. mated",
       colour = "Male refractory period") +
  facet_wrap(sex_expressed~v_cat, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        axis.text = element_text(size = 12),
        legend.position = "bottom",
         panel.spacing.x = unit(6, "mm"))

P1 + P2 + plot_layout(guides = "collect") & theme(legend.position = 'bottom')

Evolutionary dynamics

Autosomes

A_data <-
  autosome_results %>%
  filter(chromosome == "A" & dominance == 1) %>% 
  group_by(parameter_space_ID) %>% 
  slice_tail() %>% 
  ungroup() %>% 
  mutate(D_prop = D*-1,
         Fate = case_when(pop_size < 10 ~ "Extinction",
                          I_prop  >= 0.9 ~ "Invades",
                          I_prop  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.008 ~ "Low",
                           .default = "High"))
  

A_heatmap <-
  A_data %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 1) + 
  geom_vline(data = A_data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 1/(1 + 0.5*1)), linetype = 2, colour = "white", linewidth = .8) +
  geom_vline(data = A_data %>% filter(sex_expressed == "Females"),
             aes(xintercept = (0.5*1)/(1 + 0.5*1)), linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(#pnw_palette("Shuksan2", n = 5)[3],
                               pnw_palette("Shuksan2", n = 5)[4],
                               pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele",
       title = "Autosomal inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

A_heatmap

X chromosome

X_data <-
  X_results %>%
  filter(chromosome == "X" & dominance == 1) %>% 
  group_by(parameter_space_ID) %>% 
  slice_tail() %>% 
  ungroup() %>% 
  mutate(D_prop = D*-1,
         Fate = case_when(pop_size < 10 ~ "Extinction",
                          I_prop  >= 0.9 ~ "Invades",
                          I_prop  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.006 ~ "Low",
                           .default = "High"))

X_heatmap <-
  X_data %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 1) + 
  geom_vline(data = X_data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 1/(1 + 0.5)), linetype = 2, colour = "white", linewidth = .8) +
  geom_vline(data = X_data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 1/(1 + 0.5*2)), linetype = 2, colour = "white", linewidth = .8) +
  geom_vline(data = X_data %>% filter(sex_expressed == "Females"),
             aes(xintercept = (0.5)/(1 + 0.5)), linetype = 2, colour = "white", linewidth = .8) +
  geom_vline(data = X_data %>% filter(sex_expressed == "Females"),
             aes(xintercept = (0.5)/(2 + 0.5)), linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(#pnw_palette("Shuksan2", n = 5)[3],
                               pnw_palette("Shuksan2", n = 5)[4],
                               pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele",
       title = "X-linked inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

X_heatmap

Z chromosome

Z_data <-
  Z_results %>%
  filter(chromosome == "Z" & dominance == 1) %>% 
  group_by(parameter_space_ID) %>% 
  slice_tail() %>% 
  ungroup() %>% 
  mutate(D_prop = D*-1,
         Fate = case_when(pop_size < 10 ~ "Extinction",
                          I_prop  >= 0.9 ~ "Invades",
                          I_prop  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.006 ~ "Low",
                           .default = "High"))

Z_heatmap <-
  Z_data %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 1) + 
  geom_vline(data = Z_data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 1/(1 + 0.5)), linetype = 2, colour = "white", linewidth = .8) +
  geom_vline(data = Z_data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 2/(2 + 0.5)), linetype = 2, colour = "white", linewidth = .8) +
  geom_vline(data = Z_data %>% filter(sex_expressed == "Females"),
             aes(xintercept = (0.5)/(1 + 0.5)), linetype = 2, colour = "white", linewidth = .8) +
  geom_vline(data = Z_data %>% filter(sex_expressed == "Females"),
             aes(xintercept = (0.5*2)/(2*0.5 + 1)), linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(#pnw_palette("Shuksan2", n = 5)[3],
                               pnw_palette("Shuksan2", n = 5)[4],
                               pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele",
       title = "Z-linked inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

Z_heatmap

W chromosome

W_data <-
  W_results %>%
  filter(chromosome == "W" & dominance == 1) %>% 
  group_by(parameter_space_ID) %>% 
  slice_tail() %>% 
  ungroup() %>% 
  mutate(D_prop = D*-1,
         Fate = case_when(pop_size < 10 ~ "Extinction",
                          I_prop  >= 0.9 ~ "Invades",
                          I_prop  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.002 ~ "Low",
                           .default = "High"))

W_heatmap <-
  W_data %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 1) + 
  geom_vline(xintercept = 0, linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(#pnw_palette("Shuksan2", n = 5)[3],
                               pnw_palette("Shuksan2", n = 5)[4],
                               pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele",
       title = "W-linked inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 1,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

W_heatmap

Y chromosome

Y_data <-
  Y_results %>%
  filter(chromosome == "Y" & dominance == 1) %>% 
  group_by(parameter_space_ID) %>% 
  slice_tail() %>% 
  ungroup() %>% 
  mutate(D_prop = D*-1,
         Fate = case_when(pop_size < 10 ~ "Extinction",
                          I_prop  >= 0.9 ~ "Invades",
                          I_prop  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.002 ~ "Low",
                           .default = "High"))

Y_heatmap <-
  Y_data %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 1) + 
  geom_vline(xintercept = 1, linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(pnw_palette("Shuksan2", n = 5)[2],
                               pnw_palette("Shuksan2", n = 5)[4],
                               pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele",
       title = "Y-linked inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 1,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

Y_heatmap

Maternally inherited cytoplasmic chromosome

C_data <-
  C_results %>%
  filter(chromosome == "C" & dominance == 1) %>% 
  group_by(parameter_space_ID) %>% 
  slice_tail() %>% 
  ungroup() %>% 
  mutate(D_prop = D*-1,
         Fate = case_when(pop_size < 10 ~ "Extinction",
                          I_prop  >= 0.9 ~ "Invades",
                          I_prop  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.004 ~ "Low",
                           .default = "High"))

C_heatmap <-
  C_data %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 1) + 
  geom_vline(xintercept = 0, linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(#pnw_palette("Shuksan2", n = 5)[3],
                               pnw_palette("Shuksan2", n = 5)[4],
                               pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele",
       title = "Maternally inherited inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

C_heatmap

Paternally inherited cytoplasmic chromosome

P_data <-
  P_results %>%
  filter(chromosome == "P" & dominance == 1) %>% 
  group_by(parameter_space_ID) %>% 
  slice_tail() %>% 
  ungroup() %>% 
  mutate(D_prop = D*-1,
         Fate = case_when(pop_size < 10 ~ "Extinction",
                          I_prop  >= 0.9 ~ "Invades",
                          I_prop  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.004 ~ "Low",
                           .default = "High"))

P_heatmap <-
  P_data %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 1) + 
  geom_vline(xintercept = 1, linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(pnw_palette("Shuksan2", n = 5)[2],
                               pnw_palette("Shuksan2", n = 5)[4],
                               pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele",
       title = "Paternally inherited inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 7)[4], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

P_heatmap

Figure 2

#met.brewer(name="VanGogh3", n=7, type="discrete")[7]
#pnw_palette("Shuksan2", n= 5, type = "discrete")[5]

Sunset <- c("#b0cbe7", "#eba07e", "#fef7c7", "#a8554e")

A_panel <-
  A_data %>% 
  select(I_prop, 
         sex_expressed, 
         v_cat, 
         D_prop,
         refractory_period_prop_cohort_alive) %>% 
  pivot_wider(id_cols = c(D_prop, refractory_period_prop_cohort_alive, sex_expressed), 
              names_from = v_cat, values_from = I_prop) %>% 
  mutate(Result = case_when(Low >= 0.9 & High >= 0.9 ~ "Fixes with low &\nhigh encounter rate",
                            Low < 0.9 & High >= 0.9 ~ "Fixes with high\nencounter rate",
                            Low >= 0.9 & High < 0.9 ~ "Fixes with low\nencounter rate",
                            Low < 0.9 & High < 0.9 ~ "Does not fix"),
         chromosome = "Autosome") %>% bind_rows(
           
           tibble(D_prop = rep(NA, 5),
                  refractory_period_prop_cohort_alive = rep(NA, 5),
                  sex_expressed = rep("Females", 5),
                  Low = rep(NA, 5),
                  High = rep(NA, 5),
                  Result = c("Does not fix",
                             "Fixes with low\nencounter rate",
                             "Fixes with high\nencounter rate",
                             "Fixes with low &\nhigh encounter rate",
                             NA),
                  chromosome = "Autosome")
         ) %>% 
  
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Result), alpha = 1) + 
  geom_vline(data = A_data %>% 
               filter(sex_expressed == "Males") %>%
               mutate(chromosome = "Autosome"),
             aes(xintercept = 1/(1 + 0.5*1)), linetype = 2, colour = "black", linewidth = .8) +
  geom_vline(data = A_data %>% 
               filter(sex_expressed == "Females") %>% 
               mutate(chromosome = "Autosome"),
             aes(xintercept = (0.5*1)/(1 + 0.5*1)), linetype = 2, colour = "black", linewidth = .8) +
  scale_fill_manual(values = Sunset) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele") +
  facet_wrap(chromosome~sex_expressed, nrow = 1,
             labeller = label_glue('{`chromosome`} expressed in {`sex_expressed`}')) +
  coord_fixed() +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        axis.text = element_text(size =12),
        legend.title = element_markdown(),
        panel.spacing = unit(1, "lines"))

X_panel <-
  X_data %>% 
  select(I_prop, 
         sex_expressed, 
         v_cat, 
         D_prop,
         refractory_period_prop_cohort_alive) %>% 
  pivot_wider(id_cols = c(D_prop, refractory_period_prop_cohort_alive, sex_expressed), 
              names_from = v_cat, values_from = I_prop) %>% 
  mutate(Result = case_when(Low >= 0.9 & High >= 0.9 ~ "Fixes with low &\nhigh encounter rate",
                            Low < 0.9 & High >= 0.9 ~ "Fixes with high\nencounter rate",
                            Low >= 0.9 & High < 0.9 ~ "Fixes with low\nencounter rate",
                            Low < 0.9 & High < 0.9 ~ "Does not fix"),
         chromosome = "X") %>% bind_rows(
           
           tibble(D_prop = rep(NA, 5),
                  refractory_period_prop_cohort_alive = rep(NA, 5),
                  sex_expressed = rep("Females", 5),
                  Low = rep(NA, 5),
                  High = rep(NA, 5),
                  Result = c("Does not fix",
                             "Fixes with low\nencounter rate",
                             "Fixes with high\nencounter rate",
                             "Fixes with low &\nhigh encounter rate",
                             NA),
                  chromosome = "X")
         ) %>%  
  
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Result), alpha = 1) + 
  geom_vline(data = X_data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 1/(1 + 0.5)), linetype = 2, colour = "black", linewidth = .8) +
  geom_vline(data = X_data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 1/(1 + 0.5*2)), linetype = 3, colour = "black", linewidth = .8) +
  geom_vline(data = X_data %>% filter(sex_expressed == "Females"),
             aes(xintercept = (0.5)/(1 + 0.5)), linetype = 2, colour = "black", linewidth = .8) +
  geom_vline(data = X_data %>% filter(sex_expressed == "Females"),
             aes(xintercept = (0.5)/(2 + 0.5)), linetype = 3, colour = "black", linewidth = .8) +
  scale_fill_manual(values = Sunset) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele") +
  facet_wrap(chromosome~sex_expressed, nrow = 1,
             labeller = label_glue('{`chromosome`} expressed in {`sex_expressed`}')) +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  coord_fixed() +
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        axis.text = element_text(size =12),
        legend.title = element_markdown(),
        panel.spacing = unit(1, "lines"))

Z_panel <-
  Z_data %>% 
  select(I_prop, 
         sex_expressed, 
         v_cat, 
         D_prop,
         refractory_period_prop_cohort_alive) %>% 
  pivot_wider(id_cols = c(D_prop, refractory_period_prop_cohort_alive, sex_expressed), 
              names_from = v_cat, values_from = I_prop) %>% 
  mutate(Result = case_when(Low >= 0.9 & High >= 0.9 ~ "Fixes with low &\nhigh encounter rate",
                            Low < 0.9 & High >= 0.9 ~ "Fixes with high\nencounter rate",
                            Low >= 0.9 & High < 0.9 ~ "Fixes with low\nencounter rate",
                            Low < 0.9 & High < 0.9 ~ "Does not fix"),
         chromosome = "Z") %>% bind_rows(
           
           tibble(D_prop = rep(NA, 5),
                  refractory_period_prop_cohort_alive = rep(NA, 5),
                  sex_expressed = rep("Females", 5),
                  Low = rep(NA, 5),
                  High = rep(NA, 5),
                  Result = c("Does not fix",
                             "Fixes with low\nencounter rate",
                             "Fixes with high\nencounter rate",
                             "Fixes with low &\nhigh encounter rate",
                             NA),
                  chromosome = "Z")
         ) %>%  
  
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Result), alpha = 1) + 
  geom_vline(data = Z_data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 1/(1 + 0.5)), linetype = 2, colour = "black", linewidth = .8) +
  geom_vline(data = Z_data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 2/(2 + 0.5)), linetype = 3, colour = "black", linewidth = .8) +
  geom_vline(data = Z_data %>% filter(sex_expressed == "Females"),
             aes(xintercept = (0.5)/(1 + 0.5)), linetype = 2, colour = "black", linewidth = .8) +
  geom_vline(data = Z_data %>% filter(sex_expressed == "Females"),
             aes(xintercept = (0.5*2)/(2*0.5 + 1)), linetype = 3, colour = "black", linewidth = .8) +
  scale_fill_manual(values = Sunset) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele") +
  facet_wrap(chromosome~sex_expressed, nrow = 1,
             labeller = label_glue('{`chromosome`} expressed in {`sex_expressed`}')) +
  coord_fixed() +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        axis.text = element_text(size =12),
        legend.title = element_markdown(),
        panel.spacing = unit(1, "lines"))

W_panel <-
  W_data %>% 
  select(I_prop, 
         sex_expressed, 
         v_cat, 
         D_prop,
         refractory_period_prop_cohort_alive) %>% 
  pivot_wider(id_cols = c(D_prop, refractory_period_prop_cohort_alive, sex_expressed), 
              names_from = v_cat, values_from = I_prop) %>% 
  mutate(Result = case_when(Low >= 0.9 & High >= 0.9 ~ "Fixes with low &\nhigh encounter rate",
                            Low < 0.9 & High >= 0.9 ~ "Fixes with high\nencounter rate",
                            Low >= 0.9 & High < 0.9 ~ "Fixes with low\nencounter rate",
                            Low < 0.9 & High < 0.9 ~ "Does not fix"),
         chromosome = "W") %>% bind_rows(
           
           tibble(D_prop = rep(NA, 5),
                  refractory_period_prop_cohort_alive = rep(NA, 5),
                  sex_expressed = rep("Females", 5),
                  Low = rep(NA, 5),
                  High = rep(NA, 5),
                  Result = c("Does not fix",
                             "Fixes with low\nencounter rate",
                             "Fixes with high\nencounter rate",
                             "Fixes with low &\nhigh encounter rate",
                             NA),
                  chromosome = "W")
         ) %>% 
  
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Result), alpha = 1) + 
  geom_vline(data = W_data %>% 
               mutate(chromosome = "W"),
             aes(xintercept = 0), linetype = 2, colour = "black", linewidth = .8) +
  scale_fill_manual(values = Sunset, guide="none") +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele") +
  facet_wrap(chromosome~sex_expressed, nrow = 1,
             labeller = label_glue('{`chromosome`} expressed in {`sex_expressed`}')) +
  coord_fixed() +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        axis.text = element_text(size =12),
        legend.title = element_markdown(),
        panel.spacing = unit(1, "lines"))

Y_panel <-
  Y_data %>% 
  select(I_prop, 
         sex_expressed, 
         v_cat, 
         D_prop,
         refractory_period_prop_cohort_alive) %>% 
  pivot_wider(id_cols = c(D_prop, refractory_period_prop_cohort_alive, sex_expressed), 
              names_from = v_cat, values_from = I_prop) %>% 
  mutate(Result = case_when(Low >= 0.9 & High >= 0.9 ~ "Fixes with low &\nhigh encounter rate",
                            Low < 0.9 & High >= 0.9 ~ "Fixes with high\nencounter rate",
                            Low >= 0.9 & High < 0.9 ~ "Fixes with low\nencounter rate",
                            Low < 0.9 & High < 0.9 ~ "Does not fix"),
         chromosome = "Y") %>% bind_rows(
           
           tibble(D_prop = rep(NA, 5),
                  refractory_period_prop_cohort_alive = rep(NA, 5),
                  sex_expressed = rep("Males", 5),
                  Low = rep(NA, 5),
                  High = rep(NA, 5),
                  Result = c("Does not fix",
                             "Fixes with low\nencounter rate",
                             "Fixes with high\nencounter rate",
                             "Fixes with low &\nhigh encounter rate",
                             NA),
                  chromosome = "Y")
         ) %>%  
  
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Result), alpha = 1) + 
  geom_vline(data = Y_data %>% 
               mutate(chromosome = "Y"),
             aes(xintercept = 1), linetype = 2, colour = "black", linewidth = .8) +
  scale_fill_manual(values = Sunset, guide="none") +
  labs(x = ~delta~'(inbreeding depression)',
       #y = "Male refractory period\n(prob. outbred life has ended)",
       y = NULL,
       fill = "Fate of _I_ allele") +
  facet_wrap(chromosome~sex_expressed, nrow = 1,
             labeller = label_glue('{`chromosome`} expressed in {`sex_expressed`}')) +
  coord_fixed() +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1, 1)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        axis.text = element_text(size =12),
        axis.text.y = element_blank(),
        legend.title = element_markdown(),
        panel.spacing = unit(1, "lines"))

Cm_panel <-
  C_data %>% 
  select(I_prop, 
         sex_expressed, 
         v_cat, 
         D_prop,
         refractory_period_prop_cohort_alive) %>% 
  pivot_wider(id_cols = c(D_prop, refractory_period_prop_cohort_alive, sex_expressed), 
              names_from = v_cat, values_from = I_prop) %>% 
  mutate(Result = case_when(Low >= 0.9 & High >= 0.9 ~ "Fixes with low &\nhigh encounter rate",
                            Low < 0.9 & High >= 0.9 ~ "Fixes with high\nencounter rate",
                            Low >= 0.9 & High < 0.9 ~ "Fixes with low\nencounter rate",
                            Low < 0.9 & High < 0.9 ~ "Does not fix"),
         chromosome = "Cytoplasmic (m)") %>% bind_rows(
           
           tibble(D_prop = rep(NA, 5),
                  refractory_period_prop_cohort_alive = rep(NA, 5),
                  sex_expressed = rep("Females", 5),
                  Low = rep(NA, 5),
                  High = rep(NA, 5),
                  Result = c("Does not fix",
                             "Fixes with low\nencounter rate",
                             "Fixes with high\nencounter rate",
                             "Fixes with low &\nhigh encounter rate",
                             NA),
                  chromosome = "Cytoplasmic (m)")
         ) %>%  
  
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Result), alpha = 1) + 
  geom_vline(data = A_data %>% 
               filter(sex_expressed == "Males") %>%
               mutate(chromosome = "Cytoplasmic (m)"),
             aes(xintercept = 0), linetype = 2, colour = "black", linewidth = .8) +
  geom_vline(data = C_data %>% 
               filter(sex_expressed == "Females") %>% 
               mutate(chromosome = "Cytoplasmic (m)"),
             aes(xintercept = 0), linetype = 2, colour = "black", linewidth = .8) +
  scale_fill_manual(values = Sunset) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele") +
  facet_wrap(chromosome~sex_expressed, nrow = 1,
             labeller = label_glue('{`chromosome`} expressed in {`sex_expressed`}')) +
  coord_fixed() +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        axis.text = element_text(size =12),
        legend.title = element_markdown(),
        panel.spacing = unit(1, "lines"))

Cp_panel <-
  P_data %>% 
  select(I_prop, 
         sex_expressed, 
         v_cat, 
         D_prop,
         refractory_period_prop_cohort_alive) %>% 
  pivot_wider(id_cols = c(D_prop, refractory_period_prop_cohort_alive, sex_expressed), 
              names_from = v_cat, values_from = I_prop) %>% 
  mutate(Result = case_when(Low >= 0.9 & High >= 0.9 ~ "Fixes with low &\nhigh encounter rate",
                            Low < 0.9 & High >= 0.9 ~ "Fixes with high\nencounter rate",
                            Low >= 0.9 & High < 0.9 ~ "Fixes with low\nencounter rate",
                            Low < 0.9 & High < 0.9 ~ "Does not fix"),
         chromosome = "Cytoplasmic (p)") %>% bind_rows(
           
           tibble(D_prop = rep(NA, 5),
                  refractory_period_prop_cohort_alive = rep(NA, 5),
                  sex_expressed = rep("Females", 5),
                  Low = rep(NA, 5),
                  High = rep(NA, 5),
                  Result = c("Does not fix",
                             "Fixes with low\nencounter rate",
                             "Fixes with high\nencounter rate",
                             "Fixes with low &\nhigh encounter rate",
                             NA),
                  chromosome = "Cytoplasmic (p)")
         ) %>%   
  
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = Result), alpha = 1) + 
  geom_vline(data = A_data %>% 
               filter(sex_expressed == "Males") %>%
               mutate(chromosome = "Cytoplasmic (p)"),
             aes(xintercept = 1), linetype = 2, colour = "black", linewidth = .8) +
  geom_vline(data = C_data %>% 
               filter(sex_expressed == "Females") %>% 
               mutate(chromosome = "Cytoplasmic (p)"),
             aes(xintercept = 1), linetype = 2, colour = "black", linewidth = .8) +
  scale_fill_manual(values = Sunset) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Fate of _I_ allele") +
  facet_wrap(chromosome~sex_expressed, nrow = 1,
             labeller = label_glue('{`chromosome`} expressed in {`sex_expressed`}')) +
  coord_fixed() +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        axis.text = element_text(size =12),
        legend.title = element_markdown(),
        panel.spacing = unit(1, "lines"))

((W_panel | Y_panel) / Cm_panel / X_panel / A_panel / Z_panel) +
  plot_layout(guides = "collect",
              axis_titles = "collect") 

Evolutionary consequences for population fitness

Autosomes

autosome_results %>%
  filter(chromosome == "A" & time == 5 & dominance == 1) %>%   
  mutate(time = "START",
         D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.008 ~ "low",
                           .default = "high")) %>% 
  bind_rows(A_data %>% mutate(time = "END")) %>% 
  pivot_wider(names_from = time, 
              values_from = pop_size,
              id_cols = parameter_space_ID) %>% 
  mutate(pop_size_change = END / START) %>% 
  select(parameter_space_ID, pop_size_change) %>% 
  left_join(A_data) %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = pop_size_change), alpha = 1) + 
  scale_fill_gradientn(colours = moma.colors("Alkalay2"), limits = c(0.6, 1.2)) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Pop. size change",
       title = "Autosomal inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 5)[2], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

X chromosome

X_pop_heatmap <-
  X_results %>%
  filter(time == 5 & dominance == 1) %>%   
   mutate(time = "START",
         D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.006 ~ "low",
                           .default = "high")) %>%  
  bind_rows(X_data %>% mutate(time = "END")) %>% 
  pivot_wider(names_from = time, 
              values_from = pop_size,
              id_cols = parameter_space_ID) %>% 
  mutate(pop_size_change = END / START) %>% 
  select(parameter_space_ID, pop_size_change) %>% 
  left_join(X_data) %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = pop_size_change), alpha = 1) + 
  scale_fill_gradientn(colours = moma.colors("Alkalay2"), limits = c(0.6, 1.2)) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Pop. size change",
       title = "X-linked inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 5)[2], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

X_pop_heatmap

Z chromosome

Z_pop_heatmap <-
  Z_results %>%
  filter(time == 5 & dominance == 1) %>%   
  mutate(time = "START",
         D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.006 ~ "low",
                           .default = "high")) %>%  
  bind_rows(Z_data %>% mutate(time = "END")) %>% 
  pivot_wider(names_from = time, 
              values_from = pop_size,
              id_cols = parameter_space_ID) %>% 
  mutate(pop_size_change = END / START) %>% 
  select(parameter_space_ID, pop_size_change) %>% 
  left_join(Z_data) %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = pop_size_change), alpha = 1) + 
scale_fill_gradientn(colours = moma.colors("Alkalay2"), limits = c(0.6, 1.2)) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Pop. size change",
       title = "Z-linked inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 5)[2], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

Z_pop_heatmap

W chromosome

W_pop_heatmap <-
  W_results %>%
  filter(time == 5 & dominance == 1) %>%   
  mutate(time = "START",
         D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.002 ~ "low",
                           .default = "high")) %>%  
  bind_rows(W_data %>% mutate(time = "END")) %>% 
  pivot_wider(names_from = time, 
              values_from = pop_size,
              id_cols = parameter_space_ID) %>% 
  mutate(pop_size_change = END / START) %>% 
  select(parameter_space_ID, pop_size_change) %>% 
  left_join(W_data) %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = pop_size_change), alpha = 1) + 
scale_fill_gradientn(colours = moma.colors("Alkalay2"), limits = c(0.6, 1.2)) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Pop. size change",
       title = "W-linked inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 1,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 5)[2], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

W_pop_heatmap

Y chromosome

Y_pop_heatmap <-
   Y_results %>%
  filter(time == 5 & dominance == 1) %>%   
  mutate(time = "START",
         D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.002 ~ "low",
                           .default = "high")) %>%  
  bind_rows(Y_data %>% mutate(time = "END")) %>% 
  pivot_wider(names_from = time, 
              values_from = pop_size,
              id_cols = parameter_space_ID) %>% 
  mutate(pop_size_change = END / START) %>% 
  select(parameter_space_ID, pop_size_change) %>% 
  left_join(Y_data) %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = pop_size_change), alpha = 1) + 
scale_fill_gradientn(colours = moma.colors("Alkalay2"), limits = c(0.6, 1.2)) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Pop. size change",
       title = "Y-linked inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 1,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 5)[2], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

Y_pop_heatmap

Maternally inherited cytoplasmic chromosome

C_pop_heatmap <-
  C_results %>%
  filter(time == 5 & dominance == 1) %>%   
  mutate(time = "START",
         D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.002 ~ "low",
                           .default = "high")) %>%  
  bind_rows(C_data %>% mutate(time = "END")) %>% 
  pivot_wider(names_from = time, 
              values_from = pop_size,
              id_cols = parameter_space_ID) %>% 
  mutate(pop_size_change = END / START) %>% 
  select(parameter_space_ID, pop_size_change) %>% 
  left_join(C_data) %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = pop_size_change), alpha = 1) + 
scale_fill_gradientn(colours = moma.colors("Alkalay2"), limits = c(0.6, 1.2)) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Pop. size change",
       title = "Maternally inherited inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 5)[2], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

C_pop_heatmap

Paternally inherited cytoplasmic chromosome

P_pop_heatmap <-
   P_results %>%
  filter(time == 5 & dominance == 1) %>%   
  mutate(time = "START",
         D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.002 ~ "low",
                           .default = "high")) %>%  
  bind_rows(P_data %>% mutate(time = "END")) %>% 
  pivot_wider(names_from = time, 
              values_from = pop_size,
              id_cols = parameter_space_ID) %>% 
  mutate(pop_size_change = END / START) %>% 
  select(parameter_space_ID, pop_size_change) %>% 
  left_join(P_data) %>%
  ggplot(aes(x = D_prop, y = 1-refractory_period_prop_cohort_alive)) +
  geom_blank() +
  geom_tile(aes(fill = pop_size_change), alpha = 1) + 
scale_fill_gradientn(colours = moma.colors("Alkalay2"), limits = c(0.6, 1.2)) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period\n(prob. outbred life has ended)",
       fill = "Pop. size change",
       title = "Paternally inherited inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nEncounter rate: {`v_cat`}'),
             scales = "free") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = pnw_palette("Shuksan2", n = 5)[2], linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

P_pop_heatmap